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SUMMARY j (^2 ? 0 

A digital computer program coded in FORTRAN language is described that per- 
mits calc ula tion of real fluid state relations, thermodynamic properties, and 
transport properties of molecular hydrogen in any fixed ortho-para combination. 
The program is oriented toward application in numerical integration of heat- 
transfer and fluid-flow calculations, and results cover the temperature range 
from melting to dissociation for pressures up to 340 atmospheres (approx. 5000 
lb/sq in. abs). 

Properties are obtained by combinations of analytical and empirical formula- 
tions with tabulations of published data. No unpublished data are used or pre- 
sented; however, unsubstantiated extrapolations axe used to maintain the contin- 
uous range of results. Typical maximum errors with respect to published proper- 
ties are about 1 percent in density, the larger of 1 percent or 3 calories per 
gram in enthalpy, and 5 percent in specific heats. 

Any two state variables (p,T,v), or enthalpy and either pressure or specific 
volume may be specified as independent variables. Iterative solutions are used 
in calculating variables normally formulated as independent variables. Computa- 
tion speed increased with accuracy in the trial value of the dependent variable. 
Results, however, are independent of the trial values since only single-valued 
continuous formulations are used. 

The FORTRAN coding results in a total storage requirement of about 2100 
words and uses subroutines for the calculation of logarithms, exponentials, and 
square roots. Included in the total are 50 erasable words, 31 words used in re- 
turning results, and a table of constants using 215 words. The list of constants 
contains data that specify the dimensional set and the para-ortho composition, 
which can be prepared conveniently by an additional subroutine. 


INTRODUCTION 

Many analyses have indicated that molecular hydrogen has characteristics 
that make it attractive as a working fluid for nuclear -rocket propulsion applica- 



tions. Hydrogen used in such systems will probably be carried as a liquid in the 
paramodif ication because of boil-off problems typical of normal or ortho-hydrogen. 
Commercially available liquid hydrogen, however, may contain up to 5 or 10 per- 
cent unconverted ortho -hydrogen, and, in addition, many ground tests of systems 
or system components may utilize pressurized hydrogen in the normal composition. 

Many design and analysis problems in the field of nuclear-rocket propulsion, 
many of which are conveniently solved with the aid of electronic digital comput- 
ers, involve the flow and heat transfer of the working fluid. In this regard, a 
"library” of hydrogen properties that can be prepared to approximate any fixed 
para-ortho composition has been developed at the Lewis Research Center and is de- 
scribed herein. Where possible, the resulting properties are compared with the 
published experimental data. 

In anticipation of heat-transfer and fluid-flow computations (by techniques 
that may be numerical), a series of useful parameters has been selected and is 
described in detail. 

For fluid in the vapor state, saturation conditions intended for use with 
two-phase heat-transfer correlations of the Martinelli form (refs. 1 and 2) are 
also included. 

A digital computer program for achieving this end may have many criteria of 
"efficiency" including storage requirements, execution speed, accuracy of re- 
sults, and versatility. Versatility could include such factors as the ease of 
preparation, loading, calling, and updating. For use in numerical integration, 
an additional requirement is imposed, namely, that results must be unique, pre- 
cise, single-valued, and continuous over the range of interest. Thus, computed 
properties must have a precision that is unjustified by the accuracy of the re- 
sulting data. 

In order to attain the desired continuity and versatility with a minimal 
storage requirement, the accuracy of the computed properties has been compro- 
mised. The errors with respect to published hydrogen properties are presented 
herein to illustrate the validity of the techniques used, specifically, a series 
of analytical and empirical data fits that are differentiated, integrated, and 
combined with tabulations of data. 

The mechanics of the program are presented in detail, and estimations of ex- 
ecution speed and storage requirements are shown. 

No original data are presented herein, inasmuch as most of the source mate- 
rial originated at the National Bureau of Standards Cryogenic Engineering Labora- 
tory. Properties in the cryogenic range, temperatures below 100° K, have been 
compared with those shown by Roder and Goodwin (ref. 3), Goodwin, Roder, and 
Younglove (ref. 4), Goodwin, et al. (refs. 5 and 6), two progress reports (refs. 

7 and 8), and a compendium edited by Johnson (ref. 9). Transport properties are 
computed as recommended by Rogers, Zeigler, and McWilliams (ref. 10). At higher 
temperatures, the techniques of Woolley, Scott, and Brickwedde (ref. ll) are 
used. 
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ANALYSIS 


The analytical techniques used herein to formulate properties of various 
combinations of para- and ortho-hydrogen are presented in two parts: (l) the 

general relations used over the range of the library are shown, and (2) the vari- 
ous state equations used within the temperature range of interest are shown and 
discussed briefly. The definitions of symbols and the detailed algebra are dele- 
gated to appendixes A and B, respectively. The external characteristics of the 
computer program and an error analysis of the computed properties are considered 
in the RESULTS section, while the FORTRAN language coding is shown in appendixes 
C and D. 


General Relations 

Equations of state for a real fluid are commonly expressed with pressure p 
as the dependent variable and thus with temperature T and specific volume v 
as the independent variables. Two forms are 

p = cp(v,T) (la) 

p = ^ = pZRT (lb) 


where the compressibility factor Z is generally a dimensionless function of v 
and T. Naturally, the variables must have a dimensionally consistent set of 
units, but herein units will be specified only as constants are shown. 

Partial derivatives of the state parameters are used in computation of ther- 
modynamic and transport properties of the fluid and also are useful in solving 
the state equation for the normally independent temperature or specific volume by 
iterative procedures such as Newton* s method. Successive trial values are then 


v = v + (Sp/W t 


(to solve for v) 


( 2 ) 


T - T + (to solve for T) 


( 3 ) 


where Ap is the difference between the desired and computed pressures. The 
partial derivatives are obtained by differentiation of the state equation, or, if 
given in terms of the compressibility factor Z, 




(az/ST) v 

Z 


(4a) 

(4b) 
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RTZ 


(5a) 



(5b) 


Equations (4b) and (5b) are obtained by substitution of equation (lb). 

The state equations used herein are fitted to para-hydrogen properties at 
low temperatures and are largely independent of the composition at higher temper- 
atures. The difference between para- and ortho-hydrogen state relations is as- 
sumed herein to be small, so that errors in computation of poT relations for 
ortho -hydro gen compositions up to 15 percent are insignificari to errors 

in the para-hydrogen state equation itself. The composite on ant to the 

thermodynamic properties, however, and is considered. 


Enthalpy Hq and specific heat at constant pressure " q c bo* para- 
and ortho-hydrogen in the ideal state (pv = RT or Z = l) are u^wdatea as a 
function of temperature by Woolley, et al. (ref. 11 ) . The tables (ref. 11 ) are 
used directly in computation and are interpolated semilogarithmically for inter- 
mediate temperatures. Specific heat at constant volume c v q for gas in the 
ideal state is 


R R 


( 6 ) 


The enthalpy and specific heat at constant volume for the real fluid can be 
determined by integration of the equation of state along a constant temperature 
path from the ideal to the real state. Enthalpy, from reference 12 or from most 
thermodynamics texts, is given as 


cLH = c v dT + 




dv 


which is expressed in reference 11 as 


H T,p ■ %,() 
RT 


f f@) p *> - < z - « 

•'o 


(?) 


The change in specific heat at constant volume at constant temperature can also 
be determined from (refs. 11 and 12): 



( 8 ) 


4 



and 


IliP 


: v,0 


R 



> 


dp 


O) 


The equations of state used herein can he integrated in closed form, so that no 
numerical or stepwise computations are required in solution of equations (7) 
and ( 9 ) . 

The real fluid specific heat at constant pressure Cp can also be obtained 
by similar techniques. A substantial reduction in algebra is attained, however, 
through the use of c v and the partial derivatives, as shown in most thermody- 
namic texts: 



The transport properties, absolute viscosity p 0 , and thermal conductivity 
k q are computed from empirical equations shown in references 10 and 11* Again, 
the ideal-state values are given as functions of temperature, and the real fluid 
effects are determined from the state equation by using the partial derivatives 
(Sp/dT) v or the density and a correlation given by Rodgers, et al. (ref. 10). 
These equations are shown in appendix B. 


Relations for sonic velocity, or the speed of sound in the fluid, are de- 
rived in appendix B: 


c 



(id 


Rote that equation (ll) is simply a generalization of the familiar form for the 
sonic velocity of an ideal gas: 


c 



VrgRT 


(for pv = RT) 


For a wet vapor, the sonic velocity is 
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( 12 ) 



The Clapeyron equation (ref. 12) is a key step in the derivation, as shown in ap- 
pendix B. 


Equations of State 

Representation of hydrogen state relations from the melting point at 14° K 
to the onset of appreciable dissociation for any pressure up to 340 atmospheres 
is desired. In order to attain reasonable accuracy, three state equations are 
used, each for a "region" of the range of interest shown schematically in fig- 
ure 1 in terms of temperature-entropy coordinates. These three regions and the 
vapor region are assigned code numbers for use in program logic and are referred 
to herein for convenience. 

For all pressures at temperatures above 2000° K (3600° R) , hydrogen is con- 
sidered to be in the ideal state without regard to dissociation, that is, pv = RT 
or Z = 1. In this region, coded 0, the thermodynamic and transport properties 
are functions of temperature only. 

The state equation given by Woolley, et al. (ref. ll) is used for all pres- 
sures in the temperature range from 230° to 2000° K (414° to 3600° R). Thus, in 
region 2, the compressibility factor is of the form 

Z = e Bp+Cp2 (B17) 

The coefficients B and C are temperature dependent. For lack of a satisfac- 
tory substitute, the state equation is used over a larger temperature range than 
that recommended, namely, 273° to 600° K. 

Numerical consistency needed for digital computations is forced at the re- 
gional boundaries through the use of smoothing (averaging) results within 10° K 
of the boundary. Minor changes in the properties themselves are discussed in the 
section RESULTS. 

State relations for temperatures below 230° K are obtained from an empirical 
equation developed by the present author from the para-hydrogen data of refer- 
ence 5. Where the coefficients A, B, C, and D are density dependent, the pres- 
sure is fitted as 

p = A + BT + CT B + — (B23) 

T^ 

Unfortunately, the algebraic forms of the coefficients are cumbersome, and util- 
ity of the state equation may be limited to high-speed digital computer applica- 
tion. 


Unsaturated vapor is designated region -1 in figure 1. Vapor pressure and 
temperature are related by the equation for para-hydrogen from reference 7: 
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B 

log p = A + — + - + DT 


(B37) 


Saturation densities at a given temperature are determined from the inter- 
section of the vapor pressure and the liquid or gas isotars . That is, the 
pres sure -temperature relation is found from the vapor -pres sure equation (B37), 
and the liquid-gas state equation (B23) is then solved for the saturation densi- 
ties as functions of p and T. Subsequently, saturation enthalpies are evalu- 
ated at saturation density by using the liquid and gas relations. 

Some reduction in accuracy potential results from this technique, since sat- 
uration density and enthalpy are directly available in more accurate forms than 
those that result from gas or liquid computations herein. The compromise is ac- 
cepted because saturation conditions are uniquely continuous with conditions in 
the adjoining gas or liquid, so that smoothing is not required. A numerical 
evaluation of errors with respect to published data is included in the next sec- 
tion. 


RESULTS 

Three factors related to the utility of the hydrogen property program are 
discussed in this section: (l) specification of input data and output results, 

(2) estimates of storage requirements and execution speed, and (3) an analysis of 
the errors in computed results. The detail of the SUBROUTINE, in terms of the 
FORTRAN code are shown in appendix C. Also, a subroutine for conversion of units 
and preparation of constants is shown in appendix D. 


External Characteristics of SUBROUTINE STATE 

Consider first two factors that may influence the external programs which 
use the results of SUBROUTINE STATE: First, if trial values are needed, it is 

assumed that the external program supplies them. Second, the subroutine makes no 
provision for rejecting spurious input or checking the validity of results. 

Thus, spurious input may cause the return of incorrect output without warning. 

CALLING sequence . - The logic of the hydrogen property subroutine is pre- 
pared to allow for solution with any two of the state parameters (p, T, or v) 
known or with enthalpy and either p or v known. The particular pair of 
knowns is specified by the code number J in the FORTRAN language calling com- 
mand (refs. 13 and 14): 


CALL STATE (J) 


Roughly, this command means go to SUBROUTINE STATE, supply it with the value J, 
and on completion return to the next instruction. Of interest here, is the in- 


7 


terpretation of the call number J, as itemized in table I and summarized as fol- 
lows: 


CALL number , 
J 

Interpretation (known values) 

1 

Temper at ure and specific volume 

-1 

Enthalpy and specific volume 

2 

Pressure and specific volume 

3 

Temperature and pressure 

-3 

Enthalpy and pressure 

4 

Pressure and an arbitrary temperature 


The final combination is introduced in order to simplify the logic of calling 
programs in heat-transfer calculations. There is no change in the calculation of 
properties or restriction on the value of the temperature specified, but results 
are returned without modifying the storage occupied by so-called "bulk" proper- 


ties 


( e -g» * Tfiim - 


T wall + T 

2 


The specification of this storage is considered 


next. 


Common data list . - Input data and output results of the subroutine are com- 
municated to and from other programs by specifying a list of "common" locations. 
The particular list assigned for SUBROUTINE STATE is shown in table I. The sym- 
bol D represents the data assumed as an independent variable for the particular 
call, T indicates that trial values are required, X that results are returned, 
and, finally, -1 indicates that results are returned from vapor region calcula- 
tions.* Open blocks signify that the storage is neither modified nor inspected. 
The fluid region code number N FLUID is always used to initiate calc ula tions and, 
except for J = 4 calls, is returned to indicate the region in which results are 
computed. Where smoothing is used internally, the region of lower temperature 
will be indicated. Trial values of the dependent state variable are not re- 
quired if J > 0 and N FLUID = 0 , but all trial values must be supplied if 
J < 0. Trial values for saturation specific volumes v f and v g are always 
mandatory whenever vapor calculations may be used either intentionally or inad- 
vertantly. (The computations will accept input within the limits v > v •+ 

and v crit > v f > v crit/ 4 - ) g 

Storage and speed estimates . - The following estimates of storage require- 
ments and execution speed are based on FORTRAN II (ref. 14) and an IBM 7090 com- 
puter (ref. 13). In making the estimates, no provision for loading, writing, 
printing, and punching results is included. 


*■)(* 

Results are returned by "updating, " that is, by overwriting into storage. 
Thus, any previous values will be altered if results are returned. Conversely, 
if results are not returned (no X or -l), the previous values in storage will 
remain. 
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The core storage requirements of SUBROUTINE STATE are summarized in decimal 
numbers as follows : 


Requirement 

Number of 
words 

FORTRAN program 

1772 

Data list shown in table I 

31 

Temporary (erasable) storage 

50 

Constant list 

215 

Convergence specification 

1 

Unit type code 

1 

Total 

2070 


The total of about 2100 words does not include the required logarithmic , ex- 
ponential, or square-root subroutines, which may require a total of about 300 
words. Also not included in the storage requirement is the subroutine shown in 
appendix D for preparing the "constant" list. 

Estimates of execution speed have been obtained by setting up a calling rou- 
tine that sequences the subroutine through many entries and clocks the execution 
time required. The estimates shown in table II are grouped according to fluid 
region numbers, as shown in figure 1, and within each group, the independent var- 
iables used are indicated by the call letter J. Trial values are determined in- 
ternally by using N FLUID = 0 before the call or are supplied externally from a 
nonsystematic table as noted. Iterative solutions are converged to six signifi- 
cant figures. 

Execution times vary from 0.3 to 13 minutes for 10,000 cases, or from 0.002 
to 0.1 second per entry. As expected, the perfect gas computations are much 
faster than those in other regions and are independent of the dependent state 
variable. The longer times reflect more complicated formulations and iterative 
solutions; for example, with temperatures in the range of region 1 and a call of 
J = 3, additional calculations are required to determine the phase (liquid or 
gas) of the fluid. Similarly, for vapor (region -l) calculations, the state 
equations must be solved iteratively for both saturation boundaries; conse- 
quently, the time is about 0.08 second per entry. 

Where enthalpy is used as an independent variable, the execution speed is 
further reduced because of the additional iterative loop that is activated, but 
the execution times have not been evaluated. 


Errors in Computed Results 

In the following paragraphs, the deviations in computed parameters from pub- 
lished properties are discussed. Also, the properties unsupported by published 
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data are enumerated. The discussion proceeds as shown in figure 2 and covers the 
range from dissociation to melting temperatures. 


Two areas are not treated graphically as are the others, namely: 

(1) At high temperatures, a perfect gas in the ideal state is assumed. A 
figure is used only to define the temperatures above which dissociation is appre- 
ciable. 

(2) From 300° to 600° K, results are obtained from the formulations of 
Woolley, et al. (ref. 11) which were used to prepare the reference tabulations. 
Deviations are less than 0. 1 percent and therefore are not shown. 

Dissociation limits . - SUBROUTINE STATE will return results, based on a per- 
fect gas without dissociation, for temperatures up to 5000° K (9000° R). Conse- 
quently, the program should be limited by other means to within the range of 
valid properties. The area in question, which is indicated in figure 3, has been 
prepared from the calculations of reference 15 and shows the temperature at which 
given deviations occur in specific heat at constant pressure as a function of 
pressure. A 2— percent deviation is assumed as a criterion, perfect— gas results 
are valid to 1800° K (3240° R) at atmospheric pressure, or to 2200° K (3960° R) 
at a pressure level of 100 atmospheres. 

The state equation of Woolley, et al. (ref. ll) is used for temperatures up 
to 2000° K, but properties above 600° K are unsubstantiated , since the equation 
is recommended for use only to 600° K. The resulting properties are compared 
with ideal state results for lines of constant pressure as a f un ction of tempera- 
ture (fig. 4). Compressibility factor Z, enthalpy H, specific heat at constant 
pressure c p , and viscosity p are shown in figures 4(a), (b), (c), and (d), re- 
spectively. At 2000° K, deviations of up to 3 percent in compressibility factor 
are artificially smoothed, rather than up to 11 percent, as would be needed at 
600° K. For a pressure level of 340 atmospheres, the differences at 2000° and 
600 K are 1 and 3 percent in enthalpy, 0 and 1 percent in specific heat at con- 
stant pressure, and 0.5 and 3 percent in viscosity. For pressures less than the 
maximum considered herein, the effects and potential deviations are reduced, and, 
naturally, at low pressures the percentage differences between the real fluid and 
the ideal fluid are small. To reiterate, the state equation is extrapolated 
without substantiation as a means to avoid the alternative of completely artifi- 
cial smoothing. 

Transition of real fluid to cryogenic fluid . - The state equation of 
Woolley, et al. (ref . 11 ) is also used to lower recommended temperatures, namely, 
to 230° K (414° R). The state relation developed by the present author is used 
for both gaseous and liquid phases in the temperature range of region 1 14° to 

230° K. 

The deviations of various properties from the tabulations of reference 11 
are shown in figure 5 for the temperature range from 100° to 300° K. The com- 
puted properties are from regions 1 and 2 and include the effects of smoothing 
where applicable (fig. l). Deviations in percent are shown for lines of constant 
density, so that the tabulations of the reference can be used without interpola- 
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tion. The figures are shown for para-hydrogen only, hut normal (25 percent para, 
75 percent ortho) hydrogen deviations are similar in trend and magnitude. 

Pressure calculations agree with the reference data within 3/4 percent for 
densities up to 100 amagats* (about one-third of critical) and within 1- percent 

to high density (the curves are terminated at 340 atm). Enthalpy deviations are 
less than 1 percent for all densities up to 500 amagats, and specific heat at 

constant pressure agrees with the data of reference 11 to within 2~ percent. 

Specific heat at constant volume is also accurate to about 2— percent for the 
same range. 

The remainder of the properties, namely, sonic velocity, viscosity, thermal 
conductivity, and the partial derivative (dp/dT) T are unsubstantiated in the 
temperature range of 100° to 300° K. 

Cryogenic temperature range . - Deviations in the pressure calculations of 
SUBROUTINE STATE from the data of reference 7 are shown in figure 6(a) as a func- 
tion of temperature. The comparison, of which only the envelope of maximum error 
is shown, is based on integral values of density from 1 to 44 moles per cubic 
centimeter. The maximum positive deviation of 1.7 percent occurs at 22 to 24 K 
and pressures of more than 300 atmospheres (4500 lb/sq in.), whereas the maximum 
negative deviation of 0.8 percent occurs at 24^ K near the saturated— liquid 

boundary. A deviation of I— percent occurs near the saturated-liquid line at 
17° K (30.6° R) but represents an error of 0.06 atmosphere or 1 pound per square 
inch. 


Liquid-phase calculations are not constrained by the melting pressure bound- 
ary (ref. 7), so that properties will be incorrect for pressures higher than or 
temperatures lower than the melting boundary. 

Deviations in computed enthalpy from the data of reference 3 are shown in 
figure 6(b). Since the enthalpy passes through zero within the range of inter- 
est, the difference is shown in calories per gram or Btu per pound. Generally, 
the data fall within the shaded region that has a maximum deviation of about 

gi calories per gram, whereas the full-scale range of the figure is more than 

330 calories per gram. Somewhat larger deviations for the 15-atmosphere isobar 
near the critical temperature are not shown. These deviations are about -0.9, 
-2.7, 9, and 3.1 calories per gram at 36°, 35°, 34°, and 33° K, respectively. 
Otherwise, the deviations for the 15-atmosphere isobar fall within the shaded 
area. 


The deviations of specific heat at constant pressure (fig. 6(c)) show com- 
parison with data from references 7 and 12. For the most part, deviations fall 
within 4 percent. Deviations up to about 6 percent occur in the 15- and 20- 
atmosphere isobars near critical temperature. The value of specific heat at con- 
stant pressure in this range varies up to 35 Btu per pound per degree Rankine (or 


*1 amagat ~ 0.0056 lb/cu ft. 
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cal/( g ) ( °K) ) . Results on the gas side 
data, should, not ho in error 
the liquid side. 


of the vapor dome, not included in the 
by nearly as large an amount as those on 


Deviations in specific heat at constant volume in the cryogenic range are 
shown in figure 6(d). The comparison is made along lines of nearly constant den- 

cific h Iff raw data of reference 7. At low temperatures, the liquid spe- 
c fic heat at constant volume is as much as 5 percent high. The large negative 
deviations reflect the failure of the calculated values to increase rapidly near 
the critical point as the experimental values do. 

Deviations m specific heat at constant volume with respect to computed re- 
sults from reference . 8 are shown in figure 6(e) with envelopes of maximum error 
or gaseous and liquid par a -hydro gen. The value computed by SUBROUTINE STATE is 
lgher than the reference data of the National Bureau of Standards by as much as 
5 percent in the liquid phase. In the gas phase, results are within 3 percent 
except near the critical point. 

Results of all parameters near the critical point have been inspected in de- 
ail to verify that no zero or infinite values occur. Specific heat at constant 
pressure should be infinite at the critical point (ref. 16), but herein it in- 
creases to a large finite value. The author considers this a desirable compro- 
mise where results may be converged to several significant figures and only 
rarely to a mathematical identity. 


Other parameters in the cryogenic temperature range have not been directlv 
checked against experimental data; however, the viscosity and thermal conductiv- 
ity calculations are as recommended by Rogers^ et al. (ref. 10). 

Several computed properties in the vapor region reflect directly the corre- 
lation and fitting done by the National Bureau of Standards (e.g. , vapor pressure 
and saturation specific heats). Saturation density and enthalpy reflect the ac- 
curacy of the liquid or gas calculations shown in figure 6. For example, en- 
thalpy deviations of 1 and 2 calories per gram are shown in figure 6(b)- there- 
fore, the heat of vaporization should be evaluated within ±1 calorie per gram or 
about 1 percent. s * 


CONCLUDING REMARKS 

A computational technique is described that permits calculation of the prop- 
erties of various combinations of para- and ortho -hydrogen for use primarily in 
heat -transfer and fluid-flow problems. Logic is prepared to permit the use of 
British or metric units and several combinations of enthalpy, pressure, tempera- 
ture, or specific volume as the two independent variables. Properties are re- 
turned from a single entry to the computation routine, which operates at speed 
characterized by "several" entries per second. The various thermodynamic, state 
and transport properties are computed from approximate analytical expressions and 
data tabulations to cover the range from melting temperatures to limits imposed 
by appreciable dissociation for pressures up to 340 atmospheres. 
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Accuracy can be characterized by 1^ percent in state relations, about 2 cal- 

ories per gram or 1 percent in enthalpy , and 5 percent in specific heats. Near 
the critical point, hovever, deviations of up to 6 percent in specific heat at 
constant pressure and 20 percent in specific heat at constant volume are ex- 
pected. 


Levis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, February 12, 1963 
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APPENDIX A 


SYMBOLS 




British 

British 

Metric 

c 

sonic velocity 

ft/sec 

ft /sec 

m/sec 

C P 

specific heat at 
constant pressure 

ft/°R 

Btu/(lb)(°R) 

cal/(g)(°K) 

c v 

specific heat at 
constant volume 

ft/°R 

Btu/(lb)(°R) 

ca l/(g)(°K) 

g 

standard, accelera- 
tion due to 
gravity 

ft /sec^ 

ft/sec^ 

m/sec^ 

H 

enthalpy 

ft 

Btu/lb 

cal/g 

k 

thermal conduc- 
tivity 

Ib/(sec)(°R) 

Btu/( f t ) ( sec ) ( °R ) 

cal/( cm) ( sec ) (°K) 

p 

pressure 

lh/sq ft 

lb/sq in. 

atm 

R 

gas constant 

ft/°R 

cu ft/(sq in. )(°R) 

(cc)(atm)/(g) (°R) 

T 

temperature 

°R 

°R 

°K 

V 

specific volume 

cu ft/lb 

cu ft /lb 

cc/g 

X 

vapor quality 

compressibility 
factor , pv/RT 




Z 








P 

density 

lb/cu ft 

lb/cu ft 

g/cc 

P 

viscosity 

lb/(ft) (sec) 

lb/(ft)(sec) 

poises 

Subscripts: 




crit 

critical 





f saturated liquid 

fg saturated vapor minus saturated liquid 

g saturated vapor 
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p constant pressure 

ref reference 

s constant entropy 

sat saturation conditions 
T temperature 

y constant volume 

vap wet vapor 

p density 

0 ideal state 


APPENDIX B 


DETAILS OP NUMERICAL METHOD 

The analytical and empirical equations shown in this section are generally 
given in the form and dimensional system of the reference source. The conversion 
of units to a consistent set is shown in appendix D. 


Viscosity and Thermal Conductivity 

Viscosity for hydrogen In the ideal state is assumed independent of para- 
ortho composition and is computed for all temperatures from empirical equations 
given hy Woolley, Scott, and Brickvedde (ref. 11 ). 

AT 3 / 2 (T + B) 

^0 (T + C)(T + D) ' B1 ' 


where 

Pq absolute viscosity in the ideal state, poises 

T temperature, °K 

and 

A = 85.558X10" 6 poises/^ 1 / 2 
B = 650.39° K 
C = 1175.9° K 
D = 19.55° K 

Thermal conductivity in the Ideal state for temperatures below 700° K (1260° R) 

is computed from empirical equations given by reference 11, where the effect of 

para-ortho composition is introduced through the use of c_ n . 

p, u 

p 0 [A + BT + c 0 (C + HP)] 

0 m(l + E/T) ( B2 ) 


where 

k Q thermal conductivity in the ideal state, cal/(cm)(sec)(°K) 
T temperature, °K 

hQ absolute viscosity in the ideal state, poises (eq. (Bl)) 
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specific heat at constant pressure in the ideal state, cal/ (mole )(°K) 
molecular weight, g/mole 


c p,0 

m 

and 

A = 1.8341 cal/ (mole )(°K) 

B = -0.004458 cal/(mole)(°K 2 ) 

C = 1.1308 

D = 8.973X10" 4 , l/°K 
E = 3.2° K 

Ideal state thermal conductivity for temperatures above 1500 K is computed 
from correlations given by Rogers (ref. 10). The coefficients are those for a 
pressure level of* 100 atmospheres. 


k 



+ A4 + A^T 


10 


-7 


where 

k thermal conductivity, cal/(cm)(sec)(°K) 

T t emper at ur e , °K 

and 

A-^ = 3.6789x10® ( °K ) ( cal) /( cm) ( sec ) ( °K) 

A 2 = 1013.91° K 
A 3 = 5268° K 

A 4 = 4117 cal/( cm) ( sec) (°K) 

Ag = 6.982 cal/(cm) ( sec) (°K 2 ) 

Thermal conductivity for the temperature range from 700° to 1500° K is obtained 
by linear interpolation. The values for temperatures above about 500° K have 
been checked against those of reference 17 and agree within a few percent. 
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The real-fluid viscosity and thermal conductivity for region 2 are obtained 
from the ideal-state values by Enskog-type corrections. According to refer- 
ence 11, 

d = m 0 [1 + A(BpX) + B(BpX) 2 + C(BpX) 3 ] ( B3 

where 

P absolute viscosity, poises 

Pq absolute viscosity in the ideal state, poises 
and 

A = 0.175 
B = 0.7557 
C = -0.405 

k = kQ[l + A(BpX) + B(BpX) 2 + C(BpX) 3 ] , cal/(cm)(sec)(°K) (B4) 

where 

k thermal conductivity, cal/(cm)(sec)(°K) 

k Q thermal conductivity in the ideal state, cal/(cm)(sec) (°K) 
and 

A = 0.575 
B = 0.5017 
C = -0.204 

The term (BpX) in equations (B3) and (B4) is given as 

(BpX) = T^~^(Z - 1) (35) 


and an equivalent form is used in computation, namely. 



(B6) 


For cryogenic temperatures, the real-fluid corrections are computed from equa- 
tions given by Rogers, et al. (ref. 10) : 




+ Aoe ^ + A^p + A^p^ + 


(B7) 
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where 


p absolute viscosity, centipoises 

Hq absolute viscosity in the ideal state (eq. (Bl)) 

p density, g/cc 

and 

A-l = -2.515x10“® centipoise 
Ag = 3. 55 46x10” 48 centipoise 
A 3 = 400 cc/g 

A 4 = 4. 6237X10” 4 (centipoise) (cc/g) 

A 5 = -2. 6833x10”® (centipoise) (cc/g) 2 
A^ = 4.0719 (centipoises)(cc/g) 4 

O 

k =k Q + (A-l +A 2 p + A 3 p 2 +A 4 p 3 + A 5 p 4 +Agp 5 +A 7 p 6 +Agp 7 + A 9 p 8 )xl0 6 (B8) 

where 

k thermal conductivity, cal/(cm)(sec)(°K) 

kg thermal conductivity in the ideal state, from eq. (A2) 

and 

A^ = 1.84 cal/(cm)(sec)(°K) 

A 2 = 1102.6 [cal/(cm)(sec)(°K) ] [(cc/g)] 

A = 1.22648x10® [cal/(cm)(sec)(°K)] [(cc/g) 2 ] 
o 

A 4 = -1.15024x10 s [cal/(cm)(sec)(°K) ] [(cc/g) 3 ] 

A 5 = 4.95228X10 9 [cal/(cm)(sec)(°K) ] [(cc/g) 4 ] 

Ag = -1.16927X10 11 [cal/(cm) ( sec) (°K)] [(cc/g) 5 ] 

A ? = 1.56768X10 12 [cal/(cm)(sec)(°K)] [(cc/g) 6 ] 

Ag = -1. 12433X10 13 [cal/(cm)(sec)(°K) ] [(cc/g) 7 ] 

Ag = 3.36150X10 13 [cal/(cm)(sec)(°K) ] [(cc/g) 8 ] 
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Sonic Velocity 


The speed of sound in the general fluid, pv = ZRT, can he derived from the 
definition of the propagation rate of infinitesimal disturbance from refer- 
ence 12: 


c 


2 = 




(B9) 


where g is the standard acceleration due to gravity in consistent units. The 
partial derivative can he -written as 




(BIO) 


The ratio of specific heats can he expressed (from ref. 12) as 



(Bll) 


Substituting the relations (BIO) and (Bll) into the definition (B9) and extract- 
ing the square root yields 


c 



(B12) 


This relation is general for a fluid of the form pv = ZRT and reduces to 
the familiar form for the ideal state: 

c = -/rgRT 

The sonic velocity for a fluid in the vapor state where p = cp(T) may he de- 
rived as follows. Making the substitution of (BIO) into (B9) yields 



Assuming here that entropy is a function of temperature and specific volume 

gives 
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as /Ss\ /Ss\ ^ 

dV ~ \civ/rp \ST/ v dv 


For an isentropic process , 


0 



Substituting gives 


c 


2 _ 



From the Clapeyron equation for the vapor , 


and 




%g_ 

^fg 



c 


V 


T 


The sonic velocity in the vapor can then be expressed as either 



or 


c 



(B13) 


The partial derivative is the total derivative of the vapor pressure relation 
(eq. (15)) or (dp/dT) vap . This result has been derived by Hirschfelder , Curtiss, 
and Bird (ref. 16) for sonic velocity at the critical point. 


Equation of State, Region 0 

The hydrogen properties in region 0 (fig. l) are computed with the assump- 
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tion of a ( calorically) perfect gas with an equation of state as 


pv = RT (B14) 

Equation (B14) is solved for any dependent variable in closed form and therefore 
no trial values or partial derivatives are used. It is convenient to use the 
general form of the equations for specific heat at constant pressure and sonic 
velocity, however, and the derivatives are used as 



(B15) 



(B16) 


Since Z = 1, all derivatives of Z are zero, and specific heats, enthalpy, 
sonic velocity, viscosity, and thermal conductivity are functions of temperature 
only. 


Equation of State, Region 2 

For the intermediate temperature range (fig. l) the equation of state given 
by Woolley, Scott, and Brickwedde (ref. 11) is used. The correlating equation 
and a five-term series expansion are shown in the reference, but the D-|_ . . . D 4 
notation is introduced here for convenience: 


pv 

RT 


= Z = 


e pB+p 2 C 


- 1 + D^p + Dgp 2 + BgP 3 + D^p 4 


(B17) 


where 


D-j^ = B 

t>2 

D 2 " If + C 


•pO 

D 3 = V + BC 


B 


B 2 C 


D 4 = 24 + 2 


+ 


and 


where 


B = B-^ -1 / 4 + BgT' 3 / 4 + B 3 T" 5 / 4 
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Bi = 0.0055478 (degrees Kelvin^A/amagat 
Bg = -0.036877 (degrees Kelvin ) 3 / 4 /amagat 
B 3 = -0.22004 (degrees Kelvin ) 5 A/amagat 

C = C-jT" 3 / 2 + c 2 t -2 

C^_ = 0.004788 (degrees Kelvin ) 3 / 2 /amagat 2 
C2 = -0.04053 (degrees Kelvin) 2 /amagat 2 


The derivatives are readily obtained from either the exponential form or its 
series expansion. Because of the simpler algebra, the program computes partial 
derivatives from the exponential form, and because the derivatives of the expo- 
nential form cannot be integrated, the series form is used in deriving AH/RT 
and Ac v /R: 


(!) v = 


(H) T D l p + D 2 p2 + D 3 p3 + D 4 p4 

vhere 


e pB+p 2 C( B + 2 p c) |£ 

-Zp 2 (B + 2pC) (B18) 

ePB + P 2 C( pB . + p 2 c «) 

Z(pB' + p 2 C r ) 


D i 

= B’ 





D 2 

= D-jB 

' + 

c 1 



D 3 

= E 2 B 

' + 

D-^C ' 



D 4 

= b 3 b' 

1 + 

DgC' 



B' 

--fl 

fa 

1 — 1 

1 

3Bp 

+ -f- T‘ 

4 

-3/4 

C' 

II 

1 

'3^ 
, 2 

: t -3/2 

+ 2CgT" 
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The difference between the real and ideal-state enthalpy is obtained from 


h t, p ~ %,() 

RT 


f T/dZ' 
/ P\PT 

J r\ 


•) dp - (z - 1) 
h 




Then the second partial derivative of Z with respect to temperature 
constant specific volume is 


« D ”p + D^p 2 + D^p 3 + D^p 4 

\<n L 


where 


= B" 


Dg = D-jB" + D|B' + C 

= D 2 B" + DgB’ + D-jC" + C r D{ 


D'l = D 3 B" + D 3 B' + D 2 C" + C D 2 
B" = B lT -h4 ♦ § B 2 T-3A + § B 3 T-B/*) 


C ” = C 1 T “ 3 ^ 2 + 6C 2 


T“ 2 ^ 


The s 


econd integral term of equation (9) for determining Ac y /R is 
r T^/W\ dp _ f T 2 (d" + Dg P + D3P 2 + D^p 3 )dp 

l p V^Vv J 0 

- + T p2 + T p3 + T p4 ) 


The remainder of the properties are computed as shown previously. 


0 


(B20a) 

(B20b) 

at 

(B21) 


(B22) 
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Equation of State,, Region 1 


Computations in region 1 are 
pared from the para-hydrogen data 
fitting isochores with the form 


based on an empirical equation of state pre- 
of reference 5 or 7* The equation is formed by 


P = 


A + BT + CT 2 



(B23) 


The coefficients A, B, C, and D 


are shown in figures 7(a) to (d), respec- 
tively, as a function of density. The data points show the computed coefficients 
and the curves show the functions used to fit the data. The functions of densi y 
are fitted within two limitations for use in computing thermodynamic properties,, 
namely, (l) that it should he possible to integrate in closed form and to differ- 
entiate them all, and (2) that a minimum number of terms should be used to in- 
crease the confidence in the slopes. The equations used are 



where 

A-j- = -0.24337 atm/(mole/ccb'^ 7 ^ 

A 2 = 5.591X10" 10 atm/(mole/cc) 7,256 

a n = -10.67251, -0.06286125, -0.226, 0.0754 (atm) (cc) /mole 
b n = 89.507, 5.654, 25.00, 20.00 (cc/mole) 2 
p n = 16.822, 35.65, 20.00, 18.00 cc/mole 

B = pR[l + p[B x + P 1 * 5 (B 2 + B x )]j 

where 

B x = B 3 p(34.5 - p) 3 (p < 34.5 mole/cc) 

B x = B 4 (p - 34.5) 1 * 4 (p > 34.5) 

B^ = 0.027 919 58 cc/mole 

B 2 = 0.000 166 83 (cc/mole) 2,5 

B 3 = 1x10" 6 /(30x 4.5 5 ) (cc/mole) 6,5 


(B24) 


(B25) 
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B 4 = 64.5xl0- 6 /7.5 1 - 40 (cc/mole ) 3,90 


C = C lP 4 |l - C 4 p 1 / 2 (p - C 2 )|p - 34.5|J 

where 

Ci = -Ixl0~ 6 /l24 atm/(°K ) 2 (cc/mole ) 4 
Cg = 25.5 cc/mole 
C 3 = 34.5 cc/mole 

°4 = l/(l7 . 6 X -y/ 43. 1 x 8 . 6 ) (cc/mole ) 2,3 

5 


(B26) 


D = P 2 


z 


a. 


n 


Ur'tn + (P n “ ^ 

n=*l u 


J/2 


where 


a n = -8917.152, 10296.158, -371.072, 8.623, 91.596 
b n = 63.604, 76.803, 39.310, 3.684, 16.827 
P n = 5.40, 18.00, 31.30, 35.70, 39.87 


(B27) 


(atm) (cc)/ (mole) (°K 2 ) 
(cc/mole ) 2 


cc 


/mole 


The partial derivative of pressure with respect to temperature is 


T/. ■ B + 2T ° - ^ 


and the derivative with respect to specific volume becomes 


(B28) 


(B29) 


Sa 


57 = -P X< 1.970 A 3 P 1 - 970 + 7.256 A 2 p 7 * 256 


+ 2p* 


y' r n 2,3/2 - 3p3 

+ (Pn " p) J 

n— 1 u -I 


Z a n^ p n ~ p) 

[ b n + (Pn - P ) 2 ] 5/2 


> (B30) 
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SB 

37 


-p 2 R 


1 + p 


jl + p[j 


2B-i + p 1 * 5 ^^ B 


B y>]}) 


(B31) 


where 


By = B 3 [34.5 (34.5 - p) 3 - 3p (34.5 - p) 2 ] 

By = B 4 [3.5 (3p - 34. 5) 1,4 + 1.4 (p - 34. 5) 0,4 ] 


(p < 34.5 cc/mole) 
(p > 34.5) 


Sv 


-P 


4C - p 4 C 1 (c 4 p 1 / 2 )(p 


Cp) P - 34.5 



54. 5| 
34.5) 


+ p 1 / 2 


p[c 


P 4 * 5 CtC 


1°4 


P - 


. /P - C 2 

34.5 + p 


\ + P(P - C 2) ^ 


34.5 

34.5) 


( B32 


SD 


37 = 



(B33) 


Enthalpy and specific heat at constant volume are found from the ideal state 
values and the following real-gas corrections: 

7 _ PI _ _P_ = JL + JL + + D 

RT pRT pRT pR pR pRT 3 

(bZ\ _ A + _C_ 3D 
IH pRT 2 pR pRT 4 

/d 2 z\ 2A + 12D 
\3 T 2 / v pRT 5 pRT 5 


Substituting these corrections into equation (7), which was 


%,p ~ H T,0 
RT 



dp - (Z - 1) 


yields 
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Then, 



substituting into equation (9) results in 



The integrated terms, tvo of vhich are common, are obtained in closed 
equations (B24), (B26), and (B27) as 






n 


(^n + ^n 





form from 


(B34) 

< 34.5) 

> 34.5) 

(B35) 

(B36) 
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Terms in the preceding integrals and differentials are grouped to illustrate 
the similarity of the calculations in that the additional computations required 
to yield enthalpy and specific heats are less than might be expected. 


Wet Vapor Properties, Region -1 
The vapor- pressure relation given in reference 7 is 

log p = (a + ^-§-£ + DT) (B37) 

where 

p pressure, atm 
T temperature , °K 
and 


A = 2.000 620 
B = -5.009 708X10 1 °K 


C = 1.004 4 °K 
D = 1.748 495X10" 2 l/°K 


Solution of the vapor-pressure equation for unknown temperature is accom- 
plished by an iterative procedure using the slope 


dp 

dT 


= P 


~B 


(T + C)' 


+ D 


(B38) 


A trial value of temperature must be supplied. 

When the vapor pressure and temperature are known, the fluid properties at 
the saturation conditions are found by computation in region 1 on the liquid and 
then on the gas side of the vapor dome. These calculations are performed with p 
and T as independent variables and v as the dependent variable. Again, trial 
values of v_p and v are mandatory. 

Vapor quality x is then found from 


x 



(B39) 


and the following properties are assumed to vary linearly with quality as 
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H = H f + xH fg 


h = 


Hf + ^fg 


k = + xk fg 


(B40) 


Specific heat at constant volume c y for the wet vapor is derived from con- 
ditions at the saturated liquid boundary (see sketch). 



The energy change dQ is then 

dli - (sf) v aT + (57),, 

Dividing by dT and evaluating at saturation conditions yields 



By definition of specific heats, 



By virtue of the Clapeyron relations. 



Substituting and changing the variable from v to p yields the specific heat 
at constant volume for the vapor at saturated liquid density: 
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vhere c g t is computed from relations given in reference 7. The saturation- 
density equation of the same reference is differentiated with respect to tempera- 
ture; however, the correction terms recommended are omitted. From reference 7, 
then 

c + = AT(T P - T) -0 * 1 +B + CT +DT 2 +ET 3 + FT 4 + GT 5 cal/(g) (mole) (°K) 

s ax o 

(B42) 


where 

T temperature, °K 

T c critical temperature, °K 

and 

A = 1.681 574 2 cal/(g) (mole) (°K) /°K°- 9 
B = -3.280 278 9 xlO 1 cal/( g) (mole ) (°K) 

C = 6.816 987 1 cal/(g)(mole)(°K)/°K 
D = -7.319 434 1 xlO -1 cal/(g)(mole) (°K)/°K 2 
E = 3.357 435 7 XlO -2 cal/( g) (mole ) (°K) /°K 3 
F = -7.682 974 xlO -4 cal/(g) (mole ) (°K) /°K 4 
G = 6.902 922 4 xlO -6 cal/(g)(mole) (°K)/°K 5 

The saturation density, used only for the slope dp sa ^-/dT in equation (B41) is 
also obtained from reference 7: 

Psat = P c + A ( T c “ T)0 ' 38 + B ( T c " T) + C(T c - T)4/3 

+ D(T C - T) 5 / 5 + E(T C - T) 2 mole/cc (B43) 

where 

p c = 0.0155910.00005 cc/mole (ref. 7) 

T c = 32.97610.015 °K (ref. l) 

A = 0.732 346 03 XlO -2 mole/( cc) (°K) 0, 38 
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B = -0.440 742 61 xlO' 3 mole/(cc)(°K) 

C = 0.662 079 46 XlO -3 mole/( cc ) (°K)^/^ 

D = 0.292 263 63 xl0“ 3 mole/(cc)(°K) 5 / 3 
E = 0.400 849 07 XlO -4 mole/( cc) (°K) 2 

Specific heat at constant volume for the vapor at densities less than saturated 
liquid are obtained from the relation 



Isothermal integration within the vapor where the partial derivative becomes the 
total second derivative of the vapor pressure relation, a constant, yields 


c v,v c v,v=v sat 



(B44) 


Specific heat at constant pressure c p is undefined for a vet vapor, but 
the storage allocated is filled with c p of the gas at saturated vapor-density 
to avoid possible logical problems in external programs. 
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APPENDIX C 


DETAILS OP SUBROUTINE STATE 

The FORTRAN coding of SUBROUTINE STATE is shown in table III. The assign- 
ment of common storage is in some instances unique to the external programs used 
by the author, and it is anticipated that some changes would be made in applying 
the subroutine in other programs. The portions of these assignments that influ- 
ence the internal subroutine are summarized by comments in the program and are 
shown in tabular form on page 9. 

The assigned storage totals almost 300 words, and the subroutine requires 
1772 words plus that needed for subroutines for computing logarithms, square 
roots, and exponentials. 

Two parts of the subroutine should be given special attention if internal 
changes are made : 

(1) Convergence logic for J f 1 (see table i) includes artificial con- 
straint terms. 

(2) Storage assignments in JUNK are overlapping and may be conflicting if 
the order of computation is modified. 

The location of specific data tables within the major group CS is fully 
stated in the loading and unit conversion SUBROUTINE STATE S in appendix D. 
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APPENDIX D 


DETAILS OF LOADING SUBROUTINE STATE S 

The preparation of tables and constants , including the conversion of units 
from those of the references to a consistent set, and the “initialization” of 
some indexes is collected into a FORTRAN SUBROUTINE STATE S shown in table IV. 

Three alternate techniques are suggested for preparing and initializing the 
list of constants: 

(1) Load the list as “data” from previously prepared cards or magnetic tape. 

(2) Load subroutine STATE S as an additional subroutine if core storage is 
available. About 1430 storage locations are required. 

(3) Load the subroutine in a different core load. (This will require SUB- 
ROUTINE CHAIN for the IBM 7090 computer (ref. 13), and about 200 storage loca- 
tions and a tape drive unit are needed. ) 

Three sets of units are considered, namely a British set in pounds, feet, 
seconds, and degrees Rankine; a second British set in pounds per square inch and 
British thermal units; and a metric set in atmospheres, grams, centimeters, sec- 
onds, and degrees Kelvin. The conversions of these combinations into the desired 
system are shown in the coding. The first British set is requested by giving 
UNITS = 0, the second by UNITS <0, and the metric set by UNITS > 0. 

The contact with the computing subroutine is established by the assignments 
within CS, which also includes some index initialization. The ortho-hydrogen 
composition is specified by the word COMP in a common location. 

SUBROUTINE STATE S Is prepared without regard to execution speed or storage 
requirements and is not intended for repetitious execution. 
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TABLE I. - USE OP COMMON STORAGE ASSIGNMENTS IN LOGIC OF SUBROUTINE STATE 
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TABLE II. - ESTIMATES OF EXECUTION SPEED ON IBM 7090 COMPUTER. 
CONVERGENCE TO SIX SIGNIFICANT FIGURES 


Resulting Independent Trial values used 
fluid variables 

region CALL(j) 

(N FLUID) 

0 V,T (l) None required 

P,V (2) 

P,T (3) 

2 V,T (l) None required 

P,T (3) N FLUID = 0 

(none required) 

Previous result 
(ordered sequence 
■with steps of 10 
percent or more) 

1 V,T (l) None required 3- 

P,T (3) N FLUID = 0 

(none required) 3- 

Previous result 
(ordered sequence 
with steps of 10 
percent or more) 3- 


Total time Average 

10,000 entries, execution 
min time. 


0.3±0.2 


1 . 0 ± 0.2 
1 . 8 + 0. 2 

1.810.2 


2.210.2 

11.9+0.2 

11.610.2 


0.002 


0.006 

.01 


0.015 

.07 


Random 

3.012.0 

temperatures 3. 


Random 

13.0+2.0 

pressures 3 



Previous results are used as trial values for saturation specific 
volumes VL and VG. 
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TABLE III. - SUROUTINE STATE 


SUBROUTINE STATE (J) 

FLUID PROPERTY LIBRARY PARA-ORTHO HYDROGEN COMBINATIONS 

APPLICATION MUST BE EXTERNALLY LIMITED TO PREVENT USE WHERE 
WHERE DISCOCIATION OR FREEZING OCCURS OR ABOVE 5C0C PSI WHERE DATA 
WILL BE EXTRAPOLATED WITHOUT WARNING 

THE FOLLOWING ASSIGNMENTS ORIENT -STATE- TO THE CALLING PROGRAMS 

-CS-, -STORE-, -CONV-, AND -UNITS- REQUIRE 215, 31, I, AND 1 WORD 

RESPECTIVELY, AND ALL FOUR SHOULD BE IN -COMMON- 
-COMP-, REQUIRED IN -COMMON- BY -STATES-, IS NOT USED BY -STATE- 
-JUNK- IS EREASABLE AND NEED NOT BE IN -COMMON- 

-COM- AND -CORE- ARE INTERMEDIATE ASSIGNMENTS AND ARE NOT REQUIRED 
COMMON COM 

DIMENSION COM ( 20B25 ) » CORE(25,I3), CS(400), ST0REI50), 

1 JUNK ( 50 ) . „ , , 

FQU I VALFNCE ( CORE , COM ( 1 00 ) ) , (CS, COM (204-25) ) , ( STORE , CORE ( 5 1 )) , 

I ( JUNK,C0RE(276) ), (UNI TS,C0M(24) ) , ( CON V , COM ( 40 ) ) 

ASSIGNMENT OF INPUT - OUTPUT DATA STORAGE INTO -STORE- 
IMPLICIT ASSIGNMENTS ARE C(9), CP(10), CV(1I), (DP/DT)V(13) 

DIMENSIONAL SETS ARE NOTED IN THE LOADING ROUTINE -STATES- 

EQUIVALENCE (N FLU ID, STOKE ( 5 ) ) * ( P , S TURE ( 6 ) ) , ( T S , S TORE ( 7 ) ) , 

1 ( V , STORE ( 8 ) ) , ( H, STORE ( 12 )) , ( XQ , S TORE ( 14 ) ) , ( VL , S T OR E ( 1 5 ) ) , 

2 (VG, STORE! 16) ) , ( HL , S TORE ( 1 7 ) ) , ( HG , S TORE ( 1 H ) ) , 

3 (V F ILM,ST0RF ( 19) ) , (T F [ LM , S TORE ( 26 ) ) , ( V 1 SGOS , S TORE ( 2 7 ) ) , 

4 (VIS L , STORE ( 28 )) , (VIS G , STORE ( 29 )) , (THERM K , S TOR E ( 30 ) ) , 

5 (CP F ILM, STORE ( 31 ) ) 

GROUPINGS OF CONSTANTS — CONTACT WITH LOAD PROGRAM -STATES- 

EQU I VALENCE (GRAV.CS), ( KG » R GAS,CS(2I), { WB1 »C.S ( 3 ) ) » 

1 ( WB2 , CS ( 4 ) ) , ( WB 3 ,C S ( 5 ) ) , (WC1,CS(6)), (WC2,CS(7)), 

2 ( WVS1 ,GS( 8) ) , (WVS2,CS(9) ), (WVS3,CS(lO)), ( WVS4, C S ( 1 1 ) ) » 

3 ( WKl , C S ( 12)), (WK2,CS( 13) ) , ( WK3.CS ( 14 ) ) , ( WK4 , CS ( 1 5 ) ) , 

4 ( WK5 , G S ( 16)) , (EKG1 ,CS( 17) ) , ( E KG2 , CS ( 1 8 ) > , ( EKG3 , CS ( 19 ) ) , 

6 I L KG4 , C S ( 20 ) ) , ( L KG5 , C S ( 2 1 ) ) , ( EKG6 ,CS ( 22 ) ) » 

6 (VPA,CS(23) ), (VPB,CS(24) ) , (VPC,CS<25)), (VPD,CS(26)), 

7 (T CR »CS l 27 ) ) , (V CR,CS(28))» (DT FIT,CS(29)), (TF1T1,CS(30))» 

8 ( TFIT2,CS( 31 ) ) , ( TF I T3 , C S ( 32 ) ) , ( TF I T4 , CS ( 33 ) ) , ( I NDX ,CS ( 34 ) ) 

DIMENSION HAR 11 A ( 1 ) , HAR DB(l), HAR DR ( 1 ) , HAR DAB ( 1 ) , 

1 D L 1ST ( 31 ) 

EQUIVALENCE (HAR A1,CS(35))» (HAR A2,CS(36)), (HAR AZ0,CS(37)), 

2 (HAR Bl,CS(38) ), (HAR B2,CS(39)), (HAR B3,CS(40))» 

3 (HAR B4 , C S ( 4 1 ) ) , (HAR B1A,CS(42)), (HAR B2A,CS(43)), 

4 (VP LN , C S ( 44 ) ) , (VP CUN , C S ( 45 ) ) , 

5 (HAR C 1 , C S ( 46 ) ) , (HAR C2,CS(47)), (RHO FIT,CS(48)), 

6 (HAR C4 , C S ( 49 ) ) , (HAR C5,CS(50)), (HAR C6,CS(51)), 

7 (HAR C7,CS( 52) ) , ( HAR C SM , C S ( 5 3 ) ) , ( HAR 070 , C S ( 54 ) ) , 

8 (HAR D A , C S ( 5 5 ) ) , (HAR DK,CS(64)), (HAR DB,CS(73)), 

9 (HAR DAB , CS ( 82 ) ) 
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TABLE III. - Continued. SUBROUTINE STATE 


EQUIVALENCE ( DPDT ,CS ( 9 i )) , ( D2P UT2 , CS ( 92 ) ) , (JJ,CS(93)), 

1 (RHOS A,CS(94)>, (RHOS B.CSI95)), (RHOS C»CS(96)>, 

2 (RHOS D»CS(97))t (RHOS E,CS(98)), 

3 ( CSAT A,CS(99)), (CSAT B,CS(100)>, (CSAT C,CS(101)>, 

A (CSAT D,CS(102>), (CSAT E,CS(103)>, (CSAT H,CS(104)), 

5 (C SAT G »CS ( 105 ) ) , (ROD VS 1 , CS ( 106 ) ) , (ROD VS2 »CS ( 107 ) ) , 

6 (ROD VS3,CS( 108) ) , (ROD VS4 ,CS ( 109 ) ) , (ROD VS5 , CS ( 1 10 ) ) , 

7 (ROD VS6 , CS(lll)), (ROD K1,CS(112)>, (ROD K2,CS(113)), 

8 (ROD K3,CS(U4)), (ROD K4,CS(115)), (ROD K5,CS(U6)), 

9 (ROD K6,CS(117)>, (ROD K7,CS(118)) f (ROD K8,CS(119)) 
EQUIVALENCE (ROD K9,CS(120)) f (ROD A,CS(121)>, 

1 (ROD A2 , CS ( 12 2 ) ) » (ROD A3,CS(123)I, (ROD A4,CS(124)), 

2 (ROD A5,CS( 125) ) , ( D L I ST ,C S ( 12 1 ) ) , ( T700 , CS ( 15 1 ) ) , 

3 ( T1500,CS( 152 ) ) , ( TK 1500 , CS ( 1 53 > ) , ( TK INT,CS(154)) 

EQUIVALENCING OF WORKING REGION TO TEMPORARY CORE LOCATIONS 

DIMENSION HOLD ( 1 ) , T CORE(l) 

EQUIVALENCE (T CORE , JUNK ( 26 ) ) 

EQUIVALENCE (NF, TCORE), ( P 1 , TCORE ( 2 ) ) , ( T 1 , TCORE ( 3 ) ) , 

1 ( VI » TCORE (4) ) , (Cl, TC0REI5) ) , ( CP 1 , TCORE ( 6 ) ) , ( CV 1 , TCORE ( 7 ) ) , 

2 ( H 1 , T CORE ( 8 ) ) , ( P TV 1 , TCORE ( 9 ) ) , ( VI S, TCORE ( 10 ) ) , 

3 ( T HK , HOLD , TCORE ( 11 ) ) 

OVCRLAPI NG STORAGE ASSIGNMENTS 

EQUIVALENCE ( NFN 1 , TCORE ( 22 )) , ( NON, DVX , TCORE ( 2 3 ) ) , 

2 (RGT, TC0RE(24) ) , ( T TC , TCORE ( 25 ) ) 

EQUIVALENCE (N NEC, JUNK), ( NCN , J UNK ( 2 ) ) , ( NFN , JUNK ( 3 ) ) , 

1 ( J2 , JUNK ( 4 ) ) » ( J3,JUNK(5) ), ( RHO, JUNK ( 6 )) , 

2 (RHO SQ , DD2 , JUNK(7)), (RHO R,JUNK(0)), 

3 ( TX 1 , DP , T ERM D1,RH0 C2 » T ERM B,JUNK(9)), 

4 ( T X2 , DV , T ERM D2.TERM C.JUNK(IO)), 

5 ( 0V1 , VP ,C 1 TX2 » B2TX2, TERM D3, TERM D,CR L I MT , JUNK ( 1 1 )) , 

6 (C2TX1,B3TX1,BSQ2,04,Z1,D2D2,DADV, VFG , JUNK (12) » , 

7 ( WC, DD3, HFG, HAR A,JUNK(13)), (X,DWC,A INT,DT15, JUNK ( 14) ) , 

8 <DX,D2C,RRT,N,TERM A 1 , JUNK ( 15 ) ) , 

9 ( Y1,TX4,D2,SZ2TV, TERM B 1 , HAR B , TVC , JUNK( 16)) 

EQUIVALENCE ( WB , D2D3 , TERM B2 » DBDV , BTVC2 , JUNK ( 1 7 ) ) , 

1 (TERM Cl ,RHODR, TERM A2 , PVT , JUNK ( 1 8 ) ) , ( DD4, SZTV , HAR C,JUNK(19)), 

2 (BTVC.THK L,DCDV,DRS DT , DT , JUNK ( 20 )) , 

3 (BRX.DWB, RH04 , RH04C 1 , C I NT , TTC 3 , JUNK ( 2 1 ) ) , 

4 ( D3, D2D4, RT RHO,DCV R,C S A T , TSQ , JUNK ( 22 ) ) , 

5 ( DH , RHO C 3 , RHO RF T ,HAR D,JUNK(23)», 

6 ( RH01 5 , RHO C4,D INT , DEL VS , D2B , JUNK ( 2 4 ) ) , 

7 (ABS RC3.ABS RFT , DDDV , DEL TK,JUNK(25)) 

INITIATE CONTACT WITH CALLING PROGRAM 
JJ = J 

IF (UNITS) 1,3,1 

1 DO 2 1=6,31 

2 ST ORE ( I ) = STORE ( I ) * D LIST(I) 

3 PI = P 

NF = N FLUID 
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TABLE III. - Continued, 


SUBROUTINE STATE 


N NEG = -5 
N OUT = 100 
IF (JJ - 4) 4,5,5 

4 T1 = TS 
VI = V 
GO TO 9 

5 Tl = T FILM 
VI = V FILM 

9 NDN = 0 
7 NFN = NF 

J2 = XABSF(JJ) - 2 

10 J3 = J2 - 1 
NFN 1 = NFN - 1 
RGT = R GAS * T I 
TTC = T CR - T1 
NCN = 100 

IF I NFN ) 100,11,170 

COMPLETE CALCULATIONS FOR REGIONS OF PERFECT STATE 

11 IF ( J 2 ) 12,13,14 

12 PI = RGT / VI 
GO TO 15 

13 Tl = PI * VI / RG 
GO TO 15 

14 VI = RGT / PI 

15 PTV1 = RG / Vl 
PVT = - Pi / VI 
DH = 0. 

DC V R = 0. 

IF (NDN) 700,16,700 

16 IF { T 1 - T FIT 4) 17,700,700 

17 IF (Tl - T FIT 2) 601,602,602 

VAPOR REGION DETERMINATION OF VAPOR PRESSURE 

100 IF (TTC) 601,601,102 
102 TVC = Tl + VP C 
BT VC = VP B / TVC 

VP = 10.** (VP A + BTVC + Tl * VP D) * VP CON 
IF (NFN) 1C4, 130, 130 

104 BTVC = BTVC / TVC 

BTVC2 = VP LN * (VP D - BTVC) 

DPDT = VP * BTVC2 

D2P DT 2 = VP * (2. * VP LN * BTVC / TVC + BTVC2**2) 
IF ( J2 ) 109,105,106 

106 IF (JJ) 105,145,145 

105 DP = P - VP 

Tl = Tl + DP/DPDT 

IF ( ABSF ( DP ) /P - CONV ) 110,110,102 
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TABLE III. - Continued. SUBROUTINE STATE 


C 

C 

c 


OBTAIN LIQUID-SIDE SPECIFIC VOLUME AND CHECK VALIDITY 
RETURN TO LIQUID (+1) IF VOLUME IS LESS THAN SATURATED LIQUID 

109 P = VP 

110 J2 = 2 
VI = VL 
NF = -2 

N NEG = l 
GO TO 611 


C 

115 VL = VI 
VI = V 

IF (DV + CONV ) 601,116,116 
C 

C OBTAIN GAS-SIDE SPECIFIC VOLUME AND CHECK VALIDITY 

C RETURN TO GAS (+1) IF VOLUME IS MORE THAN SATURATED VAPOR 

C 

116 NF = -1 
VI = VG 
Cl = PVT 
GO TO 611 


C 


C 

C 

C 

C 


120 VG = VI 
VI = V 

IF (DV - CONV) 125,125,601 
125 VFG = VG - VL 

XQ * (V - VL) / VFG 
HL = HOLD ( 8 ) 

HG = HI 

HFG = HG - HL 

HI = HL ♦ HFG * XO 

IF (JJ) 750,780,780 

PROHIBIT ATTEMPTS TO USE V AS DEPENDENT VARIABLE 
USED IN RE-STARTING NON-VAPOR CALCS BELOW TCR 

130 DP * (VP - PI) / PI 

IF (ABSF(DP) - CONV) 700,700,131 

131 IF (DP * (VI - V CR ) ) 1 AO , 700 ,700 

140 IF ( J 2 ) 599,599,141 

141 IF (JJ) 598,145,145 
145 N NEG = 1 


NDN 

= - 

NF 

* 1 

VI 

= VG 

IF 

(VP 

146 VI 

= VL 

GO 

TO 7 


IN VAPOR REGION 


C BRANCHING TO REGIONS OF GENERAL FLUID P-T-V FORMULATION 

C 

170 RHO = l./Vl 

RHO SQ = RH0**2 
RHO R = RHO * RG 
IF (NFN1) 200,200,300 
C 
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TABLE III. - Continued. SUBROUTINE STATE 


REGION I 

CRYOGENIC TEMPERATURE RANGES, LIQUID ANO GAS 
STATE EQUATION OF HARRY 

200 RHO RFT = RHO - RHO FIT 
RT RHO = SORT F (RHO) 

ABS RFT = ABSF (RHO RFT) 

RHO 4 = RHO SQ **2 

TERM Cl = HAR C7 - RHO * HAR C6 + RHO SQ * HAR C5 
IF (RHO RFT) 201,201,202 

201 TERM Bi = HAR B3 * RHO * ABS RFT**2 
TERM B2 = TERM Bl * (4.5 * ABS RFT - 3. * RHO) 

GO TO 203 

202 TERM Bl = HAR B4 * ABS RFT**. 40 
TERM B2 = TERM Bl * (3.5 * ABS RFT ♦ 1.40 * RHO) 

TERM Cl = HAR CSM / RH04 - TERM Cl 

203 RHO 15 = RHO * RT RHO 

DBOV = 1. + RHO * ( HAR B1A + RHO 15 * (HAR B2A ♦ TERM B2)) 

HAR B = l . ♦ RHO * (HAR Bl * RHO 15 * ( HAR B2 + ABS RFT * TERM Bl)) 

RH04C1 = RH04 * HAR Cl 
RHO C 2 = RHO - HAR C2 
C RHO C3 = RHO - HAR C3 

C ABS RC3 = ABSF (RHO C3) 

RHO C4 = RT RHO * HAR C4 

HAR C = RH04C1 * (1. - RHO C4 * RHO C2 * ABS RC3) 

DCDV = 4. * HAR C - RH04C1 * RHO C4 * (ABS RC3 * (.5 * RHO C2 + 

1 RHO) + RHO C2 * S I GNF ( RHO , RHO C3)) 

C INT = RH04CI * VI * (.3333333 ♦ RHO C4 * TERM Cl) 

C 

HAR A = 0. 

DADV = 0. 

A INT = - HAR AZO 

HAR D = 0. 

DDOV = 0. 

D INT = - HAR DZO 
N = 9 

204 RHO DR = RHO - HAR DR ( N ) 

TERM D1 = 1. / (HAR DB (N) + RHO DR**2) 

TERM D2 = SQRTF (TERM D1 ) 

TERM D3 * HAR DA (N) * TERM D1 * TERM D2 
TERM D 1 = TERM D1 * RHO UR » TERM D3 

TERM D2 = TERM D2 * RHO DR * HAR DAB ( N ) 

IF (N - 5) 205,205,206 

205 HAR D = HAR D + TERM D3 
DDOV = DDD V + TERM D1 

D INT = 0 INT ♦ TERM D2 
GO TO 209 

206 HAR A = HAR A + TERM D3 
DADV = DADV 4- TERM D1 

A INT = A INT + TERM D2 
209 N = N - 1 

IF (N) 207,207,204 

207 DDOV = RHO SQ * (2. * HAR D - 3. * RHO * DDDV) 

HAR D = HAR D * RHO SQ 
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TABLE III. 


Continued. 


SUBROUTINE STATE 


C 

c 


TERM A 1 = HAR A 1 * RHO»*1.970 
TERM A2 = HAR A2 * RHO»*7.256 
DADV = RHO SQ * (2. * HAR A - 
1 + 7.256 * TERM A2 

HAR A = HAR A * RHO SQ ♦ TERM 
1.03093 = 1. / .970 
.15985 = 1. / 6.256 
A INT = A INT / 3. ♦ VI * (1 


3. * RHO * DADV) + 1.970 
A 1 + TERM A2 


03093 * TERM A1 ♦ .15985 


C 

208 RRT = RHO R * T1 

TERM B = HAR B * RRT 
TSQ = T1 * T1 
TERM C = HAR C * TSQ 
TERM D * HAR D / TSQ 

PI = HAR A + TERM B ♦ TERM C + TERM D 

PVT = - RHO * (DADV + DBDV * RRT + DCDV * TSQ + DDDV / 
PTV1 = (TERM B + 2. * (TERM C - TERM D)) / T1 
GO TO 305 


C 

210 TERM C = C INT * TSQ 
TERM D = D INT / TSQ 

OH = PI * VI - RGT ♦ A INT - TERM C ♦ TERM D 
DCV R = - 2. • (TERM D + TERM C) / RGT 

DEL VS * ROD VS1 + ROD VS2 * EXPF(RHO • ROD VS3) + RHO 
l + RHO * (ROD VS5 + RHO SQ * ROD VS6)) 

DEL TK = ROD KJ. + RHO * (ROD K2 + RHO * (ROD K3 ♦ RHO 

1 ♦ RHO * (ROD K5 + RHO * (ROD K6 ♦ RHO * ( R0DK7 

2 (ROD K8 + RHO * ROD K9 })))>)) 

225 IF (NON) 700,226,700 

226 IF ( T 1 - T FIT 2) 227,227,602 

227 IF (TTC) 700,700,102 


MID 1 - TEMPERATURE RANGE - - WOOLLEY-SCOTT-BRICKWEDDE 
NBS RP-1932 COEFFICIENTS — B AND C 


300 TX1 * l. / Tl 
TX2 = SQRTF (TX1) 

301 C1TX2 * WC1 » TX2 
C2TX1 = WC2 » TX1 
WC = T XI * (C1TX2 + C2TX1 ) 

DWC = - ( 1 • 5*C 1 TX2 ♦ 2. *C2TX 1 ) * TX1 
D2C * TX1 * ( 3.75*C1TX2 + 6.*C2TX1) 

C 

B2TX2 = WB2 * TX2 
B3TX1 = WB3 * TXl 
TX4 = SQRTF (TX2) 

WB = TX4 * (WBl + B2TX2 + B3TX1) 
q V4 B = — TX4 * .25 *(WB1 + 3.*B2TX2 ♦ 5.*B3TX1) 

D2B = TX4 * (.3125*WB1 + 1.3125*B2TX2 ♦ 2. 8125*B3TXl ) 
8SQ2 = .5 * WB**2 
C 

c DO = 1. (SUBSTITUTED) 

C D1 = WB (SUBSTITUTED) 

D2 = 8SQ2 + WC 

D3 = WB * ( BSQ2/3. + WC) 


* TERM A 1 
* TERM A2) 


TSQ) 


* (ROD VS4 

(ROD K4 
+ RHO * 
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TABLE III. - Continued. SUBROUTINE STATE 


04 = BSQ2**2/6. + WC * (BSQ2 + .5*WC) 

C 

Z1 = RHO * ( WB + RHO *(D2 ♦ RHO *(D3 + KHO » 04))) 

PI = RHO R * Tl * ( Z 1 + 1. ) 

PTV1 = PI * T X 1 * (1. + RHO * ( DWB + RHO * OWC ) ) 

PVT = - PI * RHO SO * (VI + WB -*■ RHO • ( WC + WO) 

305 DP = P - PI 

IF (J2) 525,500,505 

DO I = OWB (SUBSTITUTED) 

310 002 = WB * DWB + DWC 

003 = 02 * DWB WB * DWC 

004 = D3 * DWB + 02 * DWC 

SZTV = RHO * ( DWB + RHO *(.5*002 + RHO * ( . 3333333*D03 + RHO • 

1 . 25*DD4 ) ) ) 

OH = RGT * ( Z 1 - SZTV) 

D2D 1 = D2B (SUBSTITUTED) 

0202 = WB * D2B + D2C + DWB * DWB 

D2D3 = 02 * D2B ♦ WB * D2C + OWB * (002 + OwC ) 

0204 = D3 * D2B + D2 * 02C ♦ DWB * D03 + DWC * 002 

SZ2TV = RHO * ( D2B + RHO *(.5*D2D2 + RHO *(.3333333*0203 ♦ 

1 RHO * .25*0204) ) ) 

DC V R = - 2. * SZTV - SZ2TV 
ENSKOGS PRESSURE CORRECTIONS 
BRX = PTV 1 / RHO R - 1. 

DEL VS * 1. + BRX * ( EKG1 + BRX *(EKG2 ♦ BRX * EKG3 ) ) 

DEL TK = 1. + BRX *(EKG4 + BRX *(EKG5 ♦ BRX * EKG6 ) ) 

VERIFY CHOICE OF FLUID REGION 

320 IF (NON) 700,321,700 

321 IF ( T 1 - T FIT 4) 322,600,600 

322 IF ( T 1 - T FIT 2) 601,700,700 

CHECK CONVERGENCE USE VARIABLE CONVERGENCE FACTOR -CONV- TO 

MAKE ADJUSTMENTS TO TRIAL VALUES 

500 T 1 = MAX1F ( T 1 + DP / PTV1, .25 * Tl) 

RGT = R GAS * Tl 
TTC = T CR - Tl 
GO TO 516 

505 IF (MIN1F (PI, - PVT)) 506,506,510 

506 PI = P 

VI = VI + S I GNF (.25 * VI, VI - V CR) 

N NEG = N NEG + 1 
IF (N NEG) 510,599,519 

510 CR LIMT = VI 
IF (TTC) 512,512,511 

511 CR LIMT = .25 * ABSF (VI - V CR) 

512 VI = VI - SIGNF (MIN1F (ABSF (DP / PVT), .25 * Vl, CR LIMT), DP) 

516 IF ( ABSF ( DP > /P - CONV) 525,525,519 

519 NCN = NCN - 1 
IF (NCN) 525,525,520 

520 IF (J2) 521,521,170 

521 IF ( NFN 1 ) 208,208,530 

525 IF ( NFN 1 ) 210,210,310 
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TABLE III. 


Continued. SUBROUTINE STATE 


CONSTRAINTS RELATED TO BAD TRIAL VALUES OF TEMPERATURE 

530 IF { TTC ) 170,531,531 

531 T1 = T CR 
GO TO 601 

RESTART SEQUENCING FLUID HAS CHANGED REGIONS 

598 V = VI 

599 NF = -1 
GO TO 9 

600 NF = 0 
GO TO 9 

601 NF = 1 
GO TO 9 

602 NF = 2 
GO TO 9 

610 NFN = 0 
GO TO 620 

611 NFN = 1 
GO TO 620 

612 NFN = 2 

620 DO 621 1=2,11 

621 HOLDI I ) = TCOREt I ) 

NDN = 1 

GO TO 10 

RETURN FROM SMOOTHING CALCULATIONS 

630 IF (NFNl) 631,631,633 

631 DX = ( HOLDI 3 ) - T FIT 3) / DT FIT 
GO TO 635 

633 DX = ( HOLDI 3 ) - T FIT 1) / DT FIT 

635 X = 1. - DX 

DO 636 1=2,11 

636 T CORE ( I ) = T CORE! I) * DX + HOLDI I) * X 
GO TO 735 

THERMODYNAMIC PROPERTIES AT IDEAL STATE SEMI-LOG INTERPOLATION 

TABLE LOOK-UP OPERATIONS 

700 X = LOGFITl) 

701 INDX = INDX - 1 

702 DX = X - CSI INDX-1 ) 

IF (DX) 701,705,703 

703 IF (CS I INDX) - X) 706,705,705 
706 INDX = INDX ♦ 1 

GO TO 702 

705 DX = DX / (CS(INDX) - CSIINDX-1)) 

X = 1. - OX 

Y1 = CSIINDX+19) * X + CSI INDX+20) • DX 

HI = I CS I I NDX + 39 ) * X + CSIINDX+60) * DX) * T1 ♦ DH 

CV1 = Y 1 + R GAS • I DCV R - 1.) 

IF I J 3 ) 728,728,710 
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C 

C TRANSPORT PROPERTIES VISCOSITY AND THERMAL CONDUCTIVITY 

C 

710 VIS = WVS1 * T1 * SQRTF(Tl) • (Tl + WVS2 ) / ((Tl ♦ WVS3) * (T1 
1 + WVS4 ) ) 

IF (TI - T700 ) 715,715,720 

715 THK = VIS * (WK1 + HK2 * Tl + Y1 *(WK3 ♦ WK4 * Tl>)/(1.+ WK5 / Tl) 
GO TO 725 

720 DT 15 = T l 500 - Tl 

IF (DT 15) 723,721,721 

721 THK = TK1500 - TK INT * DT 15 
GO TO 725 

723 THK = ROD A * EXPF ((Tl - ROD A3)**2 * ROD A2) + ROD A4 + ROD A5 

1 * Tl 

725 IF ( NFN1 ) 728,727,726 

726 VIS = VIS * DEL VS 
THK = THK * DEL TK 
GO TO 728 

727 THK = THK + DEL TK 
VIS = VIS + DEL VS 

DV = (V - VI) / MAX IF (V, .1) 

728 CPI = CVl - Tl * P I V 1**2 / PVT 

Cl = VI * SQRTF ( -CPl/CVl * PVT * GRAV ) 

IF (NF + 1) 115,120,730 
C 

C BOUNDRY REGIONS - - - SMOOTHING CHECKS 

C 

730 IF ( NDN ) 735,731,630 

731 IF (NF) 735,735,732 

732 IF ( NFN1 ) 733,733,734 

733 IF (Tl - T FIT l) 735,735,612 

734 IF (Tl - T FIT 3) 735,735,610 

735 IF ( J 3 ) 736,736,740 

736 IF (JJ) 750,750,800 
C 

C FILL IN RESULTS OF FILM PROPERTIES 

C 

740 CP FILM = CPI 
V FILM = VI 
VISCOS = VIS 
THERM K = THK 
GO TO 850 
C 

C VERIFICATION OF ENTHALPY SOLUTIONS 

C 

750 DH = H - HI 

IF (A8SF(CH) / MAX IF ( ABSF ( H) , CVl) - CONV) 751,751,752 

751 IF (NF) 780,800,800 

752 N OUT = N OUT - 1 

IF (N OUT) 751,755,755 
C 

755 IF (J2) 756,800,760 

756 DT = DH / CVl 
GO TO 775 

C 

760 IF (NF) 761,770,770 
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TABLE III. 


Concluded. 


SUBROUTINE STATE 


76 L NF = I 

DH = H - HG 
VI = VG 

IF (DH) 762,762,770 
762 DH = H - HL 

IF (DH) 764, 764, 779 
764 VI = VL 

CPI = HOLD ( 6 ) 

PTVI = H0LDI9) 

PVT = HOLD ( 5 ) 

770 DT = DH / CPI 

DVX = (VI - V CFO * SIGNF (I., DT) 

V = VI + SIGNF (MIN1F ( ABSF (DT * PTVI / PVT), MAX IF (-TTC, DVX, 

1 -.25 * DVX) , .25 * VI) , DT) 

775 TS = T1 + DT 
GO TO 4 

779 HI = H 

XO = DH / HFG 

VI = VL + XQ * VFG 
NF = -1 

780 TTC3 = TTC*». 3333333 

C SAT - C SAT A * TTC**(-.l) + C SAT B / T 1 ♦ C SA T C + T 1 * 

1 (C SAT 0 + T1 » (C SAT E + T1 * (C SAT H + T1 * C SAT GDI 
DRS DT = RHOS A * TTC**(-.62) + RHOS B + TTC3 * ( RHOS C + TTC3 * 
1 RHOS D) + TTC * RHOS E 

CVl = C SAT - DPDT * DRS DT * VL**2 + D2P DT2 * (VI - VL ) 

Cl = VI * DPDT * SQRTF (GRAV / CVl) 

CVl = CVl * T 1 
PTVI = DPDT 
VIS L = HOLD( 10) 

VIS G = VIS 

VISCOS = VIS L ♦ (VIS - VIS L) « XQ 
THK L = HOLD (11) 

THERM K = THK L + (THK - THK L ) * XQ 

MOVE WORKING RESULTS TO RETURN DATA AREA 

800 DO 801 1=1,9 

801 STORE (1+4) = T CORE (I) 

RESTORE RETURN DATA TO CALLING DIMENSIONAL SET 

850 IF (UNITS) 851,855,851 

851 DO 852 1=6,31 

852 STORE ( I ) = STORE(I) / D LIST(I) 

855 RETURN 

END( 1,0, 0,0, 0,0, 0,0, 0,1, 0,0, 0,0,0) 
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TABLE IV. 


- SUBROUTINE STATE S 


SUBROUTINE STATE S 
C 

C THIS IS THE LOAD ROUTINE FOR THE FLUID LIBRARY -STATE- 

C SET-UP ROUTINE FOR PARA-ORTHO COMBINATION AND DIMENSIONAL SET 

C 

COMMON CCM 

DIMENSION COM ( 2082 5 ) * CS1400) 

EQUIVALENCE ( CS , COM { 20425 )) , ( UN I TS , COM ( 24 » ) , ( COMP , COM ( 34 ) ) 

C 

C GROUPINGS CF CONSTANTS 

C 

DIMENSION HAR DAI L ) * HAR DB(1), HAR DRIlli HAR DAB(l), 

1 D L I S T ( 3 1 ) 

EQUIVALENCE (GRAV.CS), ( RG » R GAS,CS12I), (WB1,CS(3))» 

1 ( WB2 » CS ( 4 ) ) » ( MB3 ,CS ( 5 ) ) , (WC1,CS(6)), (WC2,CS(7)), 

2 ( WVS1 ,CS( 8) ) , ( W V S2 , C S ( 9 ) ) > ( WVS3.CSI 10 ) ) , l WVS4, CS ( 1 1 )) » 

3 (MKIfCSI 12) ) , ( WK2 , C S ( 13) )» ( WK3 , CS ( 14 ) ) , ( WK4 , CS (15))* 

4 ( WK5.CSI 16) ) t (EKGl.CSI 17) ) , ( EKG2, CS ( 1 8 ) ) , ( EKG3. CS ( 19 ) ) , 

5 ( EKG4 f CS( 20 ) ) , ( EKG5 , C S ( 2 I ) ) » ( EKG6 , CS ( 22 ) ) , 

6 ( VPA.CSf 23) ) , ( VPB.CS124) ) , ( VPC , CS ( 25 ) ) , I VPD , CS ( 26 ) ) , 

7 (T CR,CS(27))» (V CR»CS(28)), (DT FIT.CS(29)) f I TF I T 1 »CS I 30 ) ) , 

8 ( TF I T2.CS ( 31 ) ) » ( TF I T3 ,CS ( 32 ) ) , ( TF I T4 , CS ( 3 3 ) ) , ( I NDX , C S ( 34 ) ) 
EQUIVALENCE (HAR AL,CS(35)), (HAR A2,CS(36)), (HAR AZO,CS(37)), 

2 (HAR B 1 f CS ( 38 ) ) t (HAR B2,CS(39)) f (HAR B3»CS(40)), 

3 (HAR B4 * CS ( 41 ) ) » (HAR B1A.CS142)), (HAR B2A,CS(43)), 

4 (VP LN, CS ( 44 ) ) , (VP C0N,CS(45)), 

5 (HAR C1,CS(46) ), ( HAR C 2 , C S ( 47 > ) , ( HAR C3 , C S ( 48 ) ) , 

6 (HAR C4 , CS ( 49 ) ) , (HAR C5,CS(50)), (HAR C6.CSI51)). 

7 (HAR C7 f CS(52))» (HAR CSM,CS(53)), (HAR DZ0,CS(54)), 

8 (HAR DA » CS ( 55 ) ) , (HAR DR,CS(64)), (HAR DB,CS(73)), 

9 (HAR DAB *CS( 82) ) 

EQUIVALENCE (RHOS A,CS(94)) f (RHOS B.CS195)), ( RHOS C»CS(96)) t 

2 (RHOS D » C S ( 97 ) ) * (RHOS E,CS(98>), 

3 (CSAT A » C S ( 99 ) ) , (CSAT B,CS(100)), (CSAT C,CS(101>), 

4 (CSAT D»CS( 102) ) , (CSAT E,CS( 103) ) > (CSAT H.CSI104)), 

5 (C SAT G tCS ( 105 ) ) > (ROD VS I *CS ( 1 06 ) ) » (ROD VS2 *CS ( 107 ) ) » 

6 (ROD VS3 »CS ( 108 ) ) * (ROD VS4 ,CS ( 109) ) , (ROD VS5 , CS ( 1 10 ) ) , 

7 (ROD VS6, CS(lll)), (ROD K I , CS ( 1 1 2 ) ) , (ROD K2,CS(113)), 

8 (ROD K3»CS( 114 ) ) , (ROD K4,CS(115)), (ROD K5,CS(116)), 

9 (ROD K6 » CS ( 117) ) » (ROD K7,CS(1I8)), (ROD K8,CS(119)) 

EQUIVALENCE (ROD K9,CS(120)), (ROD A,CS(I21)), 

1 (ROD A2,CS(122)>, (ROD A3,CS(123)), (ROD A4,CS(124)), 

2 (ROD A5»CS(I25) ) , (D L I S T , C S l 1 2 1 ) ) , ( T 700 , C S ( 1 5 1 ) ) , 

3 ( T1500.CS ( 152) ) t ( TK 1 500 . CS ( 1 53 ) ) , (TK INT,CS(154)) 

C 

C DEFINITIONS OF FLUID PROPERTY TABLE NAMES 

C 

DIMENSION T TAB(l), CP TAB(l), H TAB(l), 

I CO TAB ( 20 ) , HO TABI20) 

EQUIVALENCE (T TAB.CS ( 155 ) ) , (CP TAB ,CS ( 1 75 ) ) , (H TAB,CS(195I) 

C 

C CLEAR STORAGE AREA -CS- OF LOADING INFORMATION 

C 

DO 1 1=1,215 

1 CS ( I ) = 0. 
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TABLE IV. - Continued. SUBROUTINE STATE S 


GRAV = 32.17 
R GAS = 766.677 

TO CHANGE THE CALL - RETURN UNIT SETS MODIFY THE CORRESPONDING 
CONVERSIONS FACTORS BELOW D LIST AND STORE LOCATIONS MATCH 

DO 7 1=6,31 

7 D LIST! I ) = 1. 

IF (UNITS) 2,5,3 

ZERO OR OPEN CODE UNIT SET POUNDS-FEET-SEC-RANKINE (LITERALLY) 

NEGATIVE CODE BRITISH SET WITH PSI, BTU, FT/SEC, BTU/LB-R 

2 D L I ST ( 6 ) = 144. 

D LIST(IO) = 778.26 
0 LIST(ll) = 778.26 
D L I ST ( 12 ) = 778.26 
0 L 1ST ( 13 ) = 144. 

D LIST ( 17 ) = 778.26 
D L I ST ( 18 ) = 778.26 
D LIST (30) = 778.26 
D L I ST ( 31 ) = 778.26 
GO TO 5 


POSITIVE CODE METRIC UNITS WITH ATM, KELVIN, CC/GR, CAL/GR, 

METERS/SEC, POISES 


3 D 

LIST ( 6 ) 


2116.22 

D 

L I ST I 7 ) 


1.8 

D 

LIST18) 

= 

1. / 62.4283 

D 

LIST ( 9 ) 

= 

1. / .3048 

D 

L I ST ( 10 ) 

= 

778.26 

0 

L I ST ( ID 

- 

778.26 

0 

L I ST ( 12 ) 

= 

778.26 * 1.8 

0 

LIST( 13) 

* 

2116.22 / 1.8 

0 

LIST( 15) 

= 

1. / 62.4283 

0 

LIST ( 16) 


1. / 62.4283 

D 

L 1ST ( 17 ) 


778.26 * 1.8 

0 

LIST! 18) 

= 

778.26 * 1.8 

D 

LIST( 19) 


1. / 62.4283 

0 

LIST ( 26 ) 

= 

1.8 

0 

L I ST ( 27 ) 

= 

.0672 

0 

L I ST ( 28 ) 

= 

.0672 

D 

L I ST ( 29 ) 

= 

.0672 

D 

LIST ( 30 ) 

= 

778.26 * .0672 

0 

L I ST ( 3 1 ) 

- 

778.26 


IDEAL STATE THERMODYNAMIC FUNCTIONS 
TEMPERATURE SCALE DEGREES KELVIN 

5 T TAB ( 1 ) = 10. 

T TAB ( 2 ) = 60. 

T TAB( 3) = 80. 

T TAB ( 4 ) = 100. 

T T AB ( 5 ) = 120. 
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TABLE IV. - Continued. SUBROUTINE STATE S 


T TAB ( 6 ) = 150. 

T T AB ( 7 ) = 200. 

T T AB ( 8 ) = 250. 

T T AB ( 9 ) = 300. 

T TAB ( 10 ) = 400. 

T T AB ( 11 ) = 500. 

T TAB( 12) = 600. 

T T AB ( 13) = 700. 

T T AB ( 14 ) = 1000. 

T TAB ( 15 ) = 1500. 

T TAB ( 16 ) = 2000. 

T T AB ( 17) = 3000. 

T T AB ( 18 ) = 4000. 

T T AB ( 19 ) = 5000. 

IDEAL STATE CP CAL / MOLE - DEG K 

PARA HYDROGEN 
CP TAB(l) = 4.968 
CP TAB ( 2 ) = 5.115 
CP TAB ( 3 ) = 5.646 
CP T AB ( 4 ) = 6.455 
CP T AB ( 5 ) = 7.204 
CP TAB ( 6 ) = 7.807 
CP TAB ( 7 ) = 7.742 
CP TAB ( 8 ) * 7.380 
CP TAB ( 9 ) = 7.152 
CP TAB ( 10 ) = 7.010 
CP T AB (11) * 6.998 
CP TAB ( 12 ) = 7.010 
CP TAB ( 1 3 ) * 7.037 
CP T AB ( 14 ) = 7.219 
CP TAB ( 1 5 ) = 7.720 
CP TAB{ 16) = 8.195 
CP TAB ( 1 7 ) = 8.859 
CP TAB ( 18) = 9.342 
CP T AB ( 19 ) = 9.748 
C ORTHO HYDROGEN 

CO TAB < 1 ) = 4.968 
CO TAB (2) = 4.969 
CO TAB (3) = 5.982 
CO TAB (4) = 5.039 
CO TAB (5) = 5.170 
CO TAB (6) = 5.487 
CO TAB (7) = 6.110 
CO TAB (8) = 6.565 
CO TAB (9) = 6.809 
CO TAB (10) = 6.963 
CO TAB (11) = 6.992 
CO TAB ( 12 ) = 7.009 
CO TAB (13) = 7.036 
CO TAB (14) = 7.219 
CO TAB (15) = 7.720 
CO TAB (16) = 8.195 
CO TAB (17) = 8.859 
CO TAB (18) = 9.342 
CO TAB (19) = 9.748 
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C 

C IDEAL STATE ENTHALPY CAL / MOLE 

C PARA HYDROGEN 


H 

TAB ( 1 ) 

= 

49.6785 

H 

TAB(2> 

= 

299. 106 

H 

TAB ( 3 ) 

= 

406.015 

H 

T AB ( 4 ) 

= 

526.837 

H 

TAB( 5) 

s 

663.752 

H 

TAB( 6) 

sr 

890.605 

H 

TAB ( 7 ) 

= 

1282*70 

H 

TAB ( 8 ) 


1 660.49 

H 

TAB < 9 ) 


2 023.16 

H 

TAB < 10) 

* 

2 729. 19 

H 

TAB( 11) 

= 

3 429.24 

H 

T AB ( 12) 

- 

4 129.48 

H 

TAB( 13) 


4 831.65 

H 

TAB( 14) 

= 

6 966.23 

H 

TAB ( 15) 

= 

10 697.20 

H 

T AB { 16) 

= 

14 679.2 

H 

TAB ( 17) 

= 

23 230.9 

H 

TABl 18 > 

= 

32 345. 

H 

TAB( 19) 

= 

41 895. 


C ORTHO HYDROGEN 


HO 

TAB 

u> 


388.327 

HO 

TAB 

(2) 

= 

636.722 

HO 

TAB 

(3) 


736.179 

HO 

TAB 

(A) 

* 

836.277 

HO 

TAB 

(5) 

= 

938.227 

HO 

TAB 

(6) 

= 

1 097.78 

HO 

TAB 

(7) 

= 

1 387.90 

HO 

TAB 

(8) 

= 

1 705.80 

HO 

TAB 

(9) 

SB 

2 040.87 

HO 

TAB 

( 10) 

= 

2 731.54 

HO 

TAB 

( 11 ) 

= 

3 429.53 

HO 

TAB 

( 12) 

= 

4 129.52 

HO 

TAB 

( 13) 


4 831.66 

HO 

TAB 

( 1 A ) 

= 

6 966.23 

HO 

TAB 

( 15) 

= 

10 697.20 

HO 

TAB 

( 16) 


14 679.2 

HO 

TAB 

(17) 

= 

23 230.9 

HO 

TAB 

(18) 


32 345. 

HO 

TAB 

( 19 ) 

= 

41 895. 


C 

C SET LOOP TO PREPARE EXECUTION TABULATIONS 

C 

COMP 1=1.- COMP 
DO 11 1=1,19 

H TAB(I) = (H TAB { 1 ) * COMP1 + HO TAB ( I ) * COMP) *778.26 / 2.0157/ 
1 T TAB (I) 

CP TAB ( I ) = (CP TABU) * C0MP1 + CO TABU) * COMP) * 778.26/2.0157 
T TAB (I) = LOGF (T TABU) * 1.8) 

11 CONTINUE 
C 

C REGION 1 — LIQUID-GAS LOW TEMPERATURE REGION 

C 

ATM = 2116.22 
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TABLE IV. - Continued. SUBROUTINE STATE S 


T = 1.8 


RHO 

* 

125 836 72 




RHO 

SQ 

= RHO * RHO 




HAR 

A 1 

= - .24337 

* 

ATM 

/ RHO* 

HAR 

A2 

= 5.591 E-10 

» 

ATM 

/ RHO* 

HAR 

B 1 

= .027 919 58 

/ RHO 


HAR 

B2 

= .000 166 83 

/ 

RHO 

**2*5 

HAR 

B3 

= 1.0 E-6 / (30. * 

4. 5**3 ) 


/ RH 

HAR 

B4 

= - .000 064 5 / 

7. 5**1. 40 


/ 

HAR 

B 1 A 

= 2. * HAR 81 




HAR 

B2 A 

= 3.5 * HAR B2 




HAR 

Cl 

= - 1. / 124.0 / 1 

000 000. 

-i 

ATM / 

HAR 

C2 

= 25.5 

* 

RHO 


HAR 

C3 

= 34.5 

* 

RHO 


HAR 

C4 

= 1. / (17.6 * SQRTF (43.1) * 8.6) 

HAR 

C5 

= 1. / 5.5 




HAR 

C 6 

= (HAR C2 + HAR C3) / 4.5 



HAR 

C7 

= HAR C2 * HAR C3 

/ 3.5 



HAR 

CSM = HAR C 3 **4.5 * 

(HAR C2 / 

63. 

- HAR 

CZ - 

ATM * RHO • T **2 




HAR 

DA 

( 1 ) = -8917.152 


CZ 


HAR 

DA 

(2) = 10296.158 


CZ 


HAR 

DA 

(3) = -371.072 


CZ 


HAR 

DA 

(4) = 8.623 


CZ 


HAR 

DA 

(5) = 91.596 


CZ 


HAR 

DA 

(6) = - 10.67251 


ATM 

* RHO 

HAR 

DA 

(7) = - .06286125 


ATM 

* RHO 

HAR 

DA 

(8) = -.226 


ATM 

* RHO 

HAR 

DA 

(9) = .0754 


ATM 

* RHO 

HAR 

DR 

( 1 ) = 5.40 


RHO 


HAR 

DR 

(2) = 18.00 


RHO 


HAR 

DR 

(3) = 31.30 


RHO 


HAR 

DR 

(4) = 35.70 


RHO 


HAR 

DR 

(5) = 39.87 


RHO 


HAR 

DR 

(6) = 16.822 


RHO 


HAR 

DR 

(7) = 35.65 


RHO 


HAR 

DR 

(8) = 20.00 


RHO 


HAR 

DR 

(9) = 18.00 


RHO 


HAR 

DB 

(1) = 63.604 


RHO 

SQ 

HAR 

DB 

(2) = 76.803 


RHO 

SQ 

HAR 

DB 

(3) = 39.310 


RHO 

SQ 

HAR 

DB 

(4) = 3.684 


RHO 

SQ 

HAR 

DB 

(5) = 16.827 


RHO 

SQ 

HAR 

DB 

(6) = 89.507 


RHO 

SQ 

HAR 

DB 

(7) = 5.654 


RHO 

SQ 

HAR 

OB 

(8) = 25.00 


RHO 

SQ 

HAR 

DB 

(9) = 20.00 


RHO 

SQ 

HAR 

AZO = 0. 




HAR 

DZO = 0. 




DO 

55 

1 = 1.9 




HAR 

DAB ( I ) = 3. * HAR DA(I) / HAR 

DB ( I 

) 


RH0**3. 90 


T*»2 / RH0**4 


/ RHO**2 . 5 


/ 99.) * 8. 
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TABLE IV. 


Continued. SUBROUTINE STATE S 


IF ( I - 5 ) 51,51,52 

51 HAR DZC = HAR DZO - HAR DAB ( 1 ) * HAR DR(I) / SGRTF (HAR OB ( I ) + 

1 HAR DR ( I ) **2 ) 

GO TO 55 

52 HAR AZO = HAR AZO - HAR DAB(I) * HAR DR(I) / SGRTF (HAR DB(I) + 

1 HAR DR ( I )**2) 

55 CONTINUE 


EQUATION — FROM FIFTH PROG REPORT 


SATURATION SPECIFIC HEAT 

CHT = 778.26 / 2.0157 
C SAT A = 1.681 576 2 
C SAT B = -3.280 278 9 El 
CHT T = CHT / T 
C SAT C = 6.816 987 1 
CHT T = CHT T / T 
C SAT D = -7.319 634 1 E-l 
CHT T = CHT T / T 
C SAT E = 3.357 435 7 E-2 
CHT T = CHT T / T 
C SAT H = -7.682 974 E-4 
CHT T = CHT T / T 
C SAT G = 6.902 922 4 E-6 

SATURATION DENSITY EQUATION 
PREPARE CONSTANTS FOR USE IN 


* CHT / T** . 9 

* CHT 

* CHT T 

* CHT T 

* CHT T 

* CHT I 

* CHT T 

— FROM FIFTH PROGRESS REPORT 
DERIVITIVE 


RHOS A » .732 346 03 E-2 * .38 * RHO / T**.38 

RHOS B = -.440 742 61 E-3 * RHO / T 

RHOS C = .662 079 46 E-3 * 4./3.* RHO / T**(4./3.) 

RHOS D = -.292 263 63 E-3 * 5./3.* RHO / T**(5./3.) 

RHOS E = .400 849 07 E-4 *2. * RHO / T**2 

CONSTANTS FOR NBS RP-1932 VIRIAL EQUATION 
FROM AMAGATS AND KELVIN 


CZ = 1000. / (62.4283 

WB1 = .0055478 

WB2 = -.036877 

WB3 = -.22004 

WC1 = .004788 

WC2 = -.04053 

CONSTANTS FOR THERMAL 

WK1 = 1.8341 
WK2 = -.004458 
WK3 = 1.1308 
WK4 = .0008973 
WK5 = 3.2 


* .089888) 

* 1.8**. 25 * CZ 

* 1.8**. 75 * CZ 

* 1. 8**1. 25 * CZ 

* 1. 8**1. 5 * CZ **2 

* (1.8 * CZ ) **2 

CONDUCTIVITY COMPUTATION 

* 778.26 / 2.0157 

* 778.26/ 1.8 / 2.0157 

/ 1.8 

* 1.8 


CONSTANTS FOR IDEAL STATE VISCOSITY COMPUTATION 


WVS1 = 85.558 * .0000001 * . 0020886/ SQRTF ( 1 . 8 ) * GRAV 

WVS2 * 650.39 * 1.8 

WVS3 = 19.55 * 1.8 
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TABLE IV. - Continued. SUBROUTINE STATE S 


WVS4 * 1175.9 * 1.8 

COEFFICIENTS FOR ENSKOG PRESSURE CORRECTION OF VISCOSITY AND K 

EKG1 = .175 
EKG2 = .7557 
EKG3 = -.405 
EKG4 = .575 
EKG5 = .5017 
EKG6 = -.204 


C 


C 

c 

c 


CRYOGENIC TEMPERATURE RANGE 

RHO = 1. / 62.4283 
RHO SQ = RHO * RHO 
CK = 778.26 • .0672 
C MU = .0672 

ROD VS1 = -2.515 E-6 

ROD VS2 = 3.5546 E-18 

ROD VS3 = 400. 

ROO VS4 = 4.6237 E-4 

ROO VS5 = -2.6833 E-3 

ROD VS6 = 4.0719 

ROD K 1 = 1.84 E-6 
CKR = CK * RHO 
ROD K2 = 1102.6 E-6 
CKR = CKR * RHO 
ROD K3 = 1.22648 
CKR = CKR * RHO 
ROD K4 * -1.15024 E2 
CKR = CKR * RHO 
ROD K5 = 4.95228 E3 
CKR = CKR * RHO 
ROD K6 = -1.16927 E5 
CKR = CKR * RHO 
ROD K7 = 1.56768 E6 
CKR = CKR * RHO 
ROD K8 = -1.12433 E7 
CKR = CKR * RHO 
ROD K9 = 3.36150 E7 


TRANSPORT CORRECTIONS — LAMS 2527 


* C MU 

* C MU 

* RHO 

* C MU * RHO 

* C MU * RHO SQ 

* C MU * RHO SQ**2 

* CK 

* CKR 

* CKR 

* CKR 

* CKR 

* CKR 

* CKR 

* CKR 

* CKR 


HIGH TEMPERATURE THERMAL CONDUCTIVITY FIT FROM ROGERS 

T 700 = 700. * T 

T1500 = 1500. * T 

ROD A = 3.6789 El / 1013.91 / SQRTF ( 2. * 3.14159) * CK 

ROD A2 = - .5 * (1. / 1013.91 )**2 / T**2 

ROO A3 = 5268. * T 

ROD A4 = 4117. E-7 * CK 

ROD A5 = 6.982 E-7 * CK / T 

TK 1500 * ROD A * EXPF (IT1500 - ROD A3)**2 * ROD A2 ) + ROD A4 
1 + ROD A5 * T 1500 

VIS 700 = WVS1 * T 700 * SQRTF C T700 ) * ( T700 ♦ WVS2 ) / 

1 ( ( T700 + WVS3 ) * ( T700 + WVS4 ) ) 
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TABLE IV. - Concluded. SUBROUTINE STATE S 


TK 700 = VIS 700 * (WK1 + WK2 * T700 + CP TAB(13) * (WK3 + WK4 * 

1 T700 ) )/(!.+ WK5 / T700 ) 

TK INT = (TK 1500 - TK 700) / (T 1500 - T 700) 

VAPOR PRESSURE RELATIONSHIPS — FROM FIFTH PROGRESS REPORT 7246 

VP A = 2.000 620 
VP B = -5.009 708 * 10. 

VP C = 1.004 4 
VP D - 1.748 495 E-2 
VP LN = LOGF { 10. ) 

VP CON = ATM 

CONSTANTS DEFINING THE REGIONAL BOUNDRIES 

T CR = 32.994 * 1.8 
V CR = 65.5 / 125.8 
T FIT 1 = 220. * 1.8 
T FIT 2 = 240. * 1.8 
T FIT 3 = 1990. * 1.8 
T FIT 4 = 2010. * 1.8 
DT FIT = 20. *1.8 

INDX = 157 
C 

CZ = COM PI * 100. 

CX = COMP * 100. 

WRITE OUTPUT TAPE 6, 56, CZ, CX 

56 FORMAT (47H0 FLUID PROPERTY LIBRARY WITH A COMPOSITION OF F5.1, 

1 20H PERCENT PARA-, AND F5.1, 35H PERCENT ORTHO- HYDROGEN, 11/1/ 
261 ) 

RETURN 

C 

END 1 1 , 0,0, 0,0, 0,0, 0,0, 1,0, 0,0, 0,0) 


* T 

* T 
/ T 
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Temperature 


Dissociation 



Rank in e 
V- 3600+19 
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Temperature 



Figure 3 - Deviations in specific heat at constant pressure due to dissociation (based on equilibrium cal- 

culations of ref. 15). 
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CD 







Viscosity, li/ng “ U Specific heat at constant 

percent pressure, c-n/c-n n - 1, 



(c) Specific heat at constant pressure. 



720 1440 2160 2880 3600 43 

Temperature, °R 


(d) Viscosity* 


Figure 4* - Concluded. Real-gas properties resulting from 
extrapolation of state equation from reference 11. 
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Deviations in computed properties of para-hydrogen "below 100' 


PSi / 

Enthalpy deviation, — , cal/g 



0 18 36 54 72 90 108 126 144 162 180 

Temperature, °R 

(h) Enthalpy from data of reference 3. Pressure: p = 1, 2, 3, 5, 7, 10, 20, 30, 40, 50, 

60, 80, 100, 120, 160, 200, 300, 340 atmospheres. 

Figure 6. - Continued. Deviations in computed properties of para-hydrogen "below 100° K. 
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Figure 6. - Continued, Deviations in computed properties of para-hydrogen helow 100' 



Figure 6. - Continued. Deviations in computed properties of para-hydrogen below 100° 
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Coefficient B, atm/°K 



(b) Coefficient B. B = PR{1 + PfBj. + P 1 ' 5 ^ + B x ])} 
p < 34.5, B x = B 3 p(34.5 - p) 3 
P > 34.5, Bx = B 4 (p - 34. 5) 1 - 4 . 

Figure 7. - Continued. Coefficients of state equation for region 1. Data of 
reference 7 fitted by P = A + BT + CT^ + d/T^. 







